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A Time for 


Reckoning 


The human mind, when it sails by dead reckoning, 
without the possibility of a fresh observation, perhaps 
without the instruments necessary to take one, will 
sometimes bring up in very strange latitudes. 


James Russell Lowell (1890), “Witchcraft” [1] 


le the mind with its hundred billion neurons is the most 
intricate amalgam in the known universe, we should expect to do 
marvelous things with it. As we form electrical pathways among 
the hundred trillion synapses in the mind when performing a task, 
or even in anticipating a task, a library of pertinent mind maps 
and templates would be invaluable. This book investigates a por- 
tion of sucha library, a shelf in the mathematics section, and browsing 
there can be an enjoyable and interesting experience for us. 

This book describes techniques of computation and approxima- 
tion that may be used to rapidly and mentally calculate mathemat- 
ical quantities, including results of arithmetic operations and 
values of elementary functions. Some of these methods are very 
old, some are relatively recent, and some are new. Overall, the 
idea is to enjoy developing these capabilities. The concept of dead 
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reckoning is meant to convey the adventurous and challenging 
spirit that I perceive in this pursuit. At the very least, we will 
experience a portion of the panorama of elementary mathematics 
that mankind has had thousands of years to develop and that 
includes many calculational methods, of which we normally 
consider the ones we learned in school to be the only (or even 
best) ones. For mental calculations, they usually are not. An 
analogous situation occurred in the development of our present 
pen-and-paper methods of calculation that minimize erasures, 
supplanting the ancient sand reckoning methods in which digits 
were continuously and instantly erased and overwritten in the 
course of the calculations. 

This does not pretend to be a book for teaching arithmetic. It 
is intended for those in high school, college, or professional 
occupations who are intrigued by mathematics, regardless of 
whether they are considered “good” at it or not. Unlike some 
books on calculational shortcuts, it is not a business math book— 
you won’t find any dollar signs here. No exercises are included; the 
world is full of them. On the other hand, this book is replete with 
examples, which usually offer more personal insight than the often 
forbidding symbolic formulas. 

For those who appreciate numbers and their interrelationships, 
the ability to perform mental calculations can provide a great deal 
of personal satisfaction and recreation. Speed is without a doubt 
a major objective of these techniques, but the sheer power of the 
mind to extract, say, square roots to sixteen or more digits offsets, 
I feel, the several minutes required for someone of my limited 
capabilities to do it. Therefore, speed is not the only concern and 
some methods are included that do take some time. Some of these 
algorithms are very interesting and can be thought-provoking in 
their own right. 

Some may say that this book is a relic, given the calculating 
power available to anyone these days. | counter that the contents 
of this book are in fact both timeless and timely. Certainly the 
enjoyment of exploring mathematics and the challenge of stretch- 
ing our personal capabilities will always exist for a cognitive 
people. On the other hand, I propose that this is a very appro- 
priate time for presenting such calculational methods, for two 
curiously interrelated reasons. First, the proliferation of electronic 
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calculators and computers throughout our lives and educational 
systems are eliminating calculational techniques from our memory 
and from our children’s education, threatening to dim our already 
narrow view of this very rich field. 

Ironically, the material in this book is also timely because a 
substantial amount of work on fast calculational algorithms was 
performed in the early days of mechanical and electronic com- 
puters. Limited in speed and memory and structured as sequential 
operations (sound familiar?), these early computers demanded 
simple, fast, and accurate approximations to mathematical func- 
tions. The work in this field, particularly that for decimal 
machines, is reflected in a variety of techniques gathered and pre- 
sented in this book. Notably, more recent work in computer 
algorithms is much less significant; numeric coprocessors, 
extended-bit precision, parallel processing, and so forth have made 
the results much less applicable to humans. Newer algorithms, 
such as the Fast Fourier Transform (FFT) method of multi- 
plication, often rely on internal architectures of microprocessors 
and their siblings, arrays of available memory locations, and/or 
inherent bit-shifting operations, none of which find analogy in our 
thought processes. 

My own inability to readily approximate elementary functions 
was first brought home several years ago while | was working in 
an engineering laboratory. I required the value of sin 28° to 
quickly check a test result. Not having a calculator with me, | 
came to the shocking realization that, excluding the tedium of 
calculating terms by hand in the slowly converging Maclaurin 
power series (in radian units, no less), | was lost. Throughout my 
physics and mathematics education, trigonometric functions were 
always found from tables, slide rules, or calculators. 

The next day, as it turned out, the need arose in a meeting to 
estimate the square root of 39. Again, a ready estimate to any 
significant accuracy at all was impossible. | was appalled. 

I was determined to do something about this. As a graduate 
teaching assistant | had quietly enjoyed calculating results of 
arithmetic expressions before students could locate their calcu- 
lators, and I felt I had a fair aptitude with numbers. | had recently 
read Steven B. Smith’s excellent book, The Great Mental Calcu- 


lators [2], and was somewhat relieved to recall that these were not 
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feats that calculating prodigies were generally capable of doing 
either. Trigonometric functions were definitely not their province, 
and non-integer roots were seldom attempted. 

Let us digress here for a moment to talk about such prodigies, 
as we will occasionally encounter them in the later chapters. At 
the outset, | would like to make it clear that | am not one of these 
lightning calculators, a very talented and practiced breed. How- 
ever, it also seems appropriate in these days of high-speed pocket 
calculators to add that despite popular notions, these lightning 
calculators, which we can certainly strive to emulate, were and are 
by no means instantaneous in deriving their results. Due to the 
paucity of reliable measurements of response times, as well as the 
large deviations evident among specific types of presentations and 
calculators, | am extremely reluctant to elaborate on this and refer 
the interested reader to Smith’s detailed exposition. I would like 
to caution that descriptions of performances by lightning calcu- 
lators are poor yardsticks. They are generally written to flatter the 
calculator, they ignore such delay tactics as conversation and 
writing or repeating numbers, and they arbitrarily use phrases such 
as “in an instant,” “in a flash,” “momentarily,” and so forth, which 
are not quantified. 

With this established up front, and only to demonstrate the 
finite time required, I point out a test, better than most, reported 
in 1894 by Alfred Binet. In the test, the time required to multiply 
6241 by 3635 ranged from 21 seconds (Jacques Inaudi) to 70.5 
seconds (Ugo Zaneboni). Pericles Diamandi took 127 seconds to 
multiply 8637 by 4538. On the other hand, Johann Dase was 
credited by a knowledgeable observer in 1861 with mentally 
multiplying two eight-digit numbers (79532853 and 93758479) 
correctly in 54 seconds. It is also true that two-digit by two-digit 
multiplications are often immediately drawn, at least with very 
little adjustment, from memory. 

Equally important is the realization that the results are not 
generated by a spontaneous process. From Hofstadter [3]: 


From such people’s introspection, as well as from extensive 
research by psychologists, it has been ascertained that noth- 
ing occult takes place during the performances of lightning 
calculators, but simply that their minds race through inter- 
mediate steps with the kind of self-confidence that a natural 
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athlete has in executing a complicated motion quickly and 
gracefully. They do not reach their answers by some sort of 
instantaneous flash of enlightenment (though, subjectively, it 
may feel that way to some of them), but—like the rest of us— 
by sequential calculation . . . 


Mental calculators have usually been extremely talented indi- 
viduals who mastered their craft through sheer diligence from an 
early age. Efforts typically involved memorization of large tables, 
including two-digit endings for root extraction, logarithmic tables, 
factors of large numbers, and powers of two-digit numbers. Cou- 
pled with an innate number sense cultivated over years, the 
calculator historically provided astonishing feats of computation of 
limited types, often without an explicit algebraic understanding of 
the techniques. On the other hand, many noted theoreticians 
were considered lightning calculators, including Pascal, Gauss, 
Euler, Ampere, Wallis (who purposefully developed the ability in 
middle age), Aitken, Fermi, and Von Neumann. 

Actually, the types of calculations performed by professional 
lightning calculators are often specialized to a degree that one 
would rarely encounter them even in technical work. These 
include extracting higher roots of exact powers, finding multiple 
squares whose sum is given (for example, every integer can be 
expressed as at least one sum of four squares), calculating com- 
pound interest, and performing day—date calendar functions. 

I began, then, to survey the methods of lightning calculators as 
well as what I discovered to be a rather extensive field in 
numerical analysis, approximation of elementary functions. This 
was simply an interest of mine initially. After a couple of years, 
I decided to write an article for a mathematics journal. However, 
the article grew quite large and the scope eventually widened to 
culminate in this book instead. I am quite sure that you will not 
find another book like this one; I looked and ended up researching 
this one. 

This book is aimed at those of us who, through our interest in 
mathematics, have acquired a general appreciation for numbers 
(i.e., a number sense), but have limited desire, time, or talent to 
develop the tools of lightning calculators. While the techniques 
and strategies presented are amenable to mental calculation, they 
are useful, of course, for quick pen-and-paper calculations as well. 
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Ideally, I think, algorithms designed to be performed mentally 
should have the following characteristics: 


1. Multiplications that are limited to two-digit by three-digit 
at most 

2. Divisions that are limited to two-digit divisors at most 

3. A short-term memory requirement of no more than two 

numbers at any one time 

A very minimum number of memorized table values 

. A relatively quick result, accurate to at least three digits, 
with the capability of extending this accuracy if so desired 


wr 


This book, then, is not meant to provide a general collection 
of approximation techniques, but rather to merge selected tools of 
the mental calculator with relatively modern techniques of 
numerical analysis. To be precise, it is comprised of a set of useful 
methods and algorithms for rapidly and mentally finding products, 
quotients, roots of nonperfect powers, trigonometric values, and 
logarithmic and exponential values. Extremely detailed error anal- 
yses are omitted here to preserve the flavor of the discussion (and, 
to be honest, my taste for the subject). If you are interested in 
more detail, please consult the references cited; ] have tried to be 
meticulous. Brief derivations of formulas, however, are generally 
given as they occur, although I trust the mathematicians among 
us will forgive my lack of rigor. 

A notation mention is unavoidable here, because in this book 
certain operations involve working with groups of digits in a 
number. The notation |n represents a two-digit number string; if 
more digits exist in n, they are “melded,” or added, to the digits 
to the left of the “|” sign. For example, 31129 = 4|29 = 429. We 
simply want to work with hundreds groups in the number, and 
these can carry or borrow as needed from neighboring groups. 

Don’t be confused by the fact that there are three digits in the 
rightmost hundreds group of 31129. This can occur in an inter- 
mediate result of calculating with individual groups; the melding 
process returns each hundreds group to two digits. To illustrate, 
consider the numbers 184 and 245. To add these, we can split 
them up into two-digit groupings (1184 and 2145) and add each 
piece separately (1184 + 2145 = 31129). We then meld the result 
into the final answer: 
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+ 
429 


where the | is a carry into the group on its left, as in addition. 

While this notation may seem unnecessary at this point, we will 
find that in multiplication, division, and square root extraction, 
we can work comfortably with separate two-digit groups within 
numbers and then meld the final answer back to a number with 
two-digit groupings, whence we drop the vertical bars. 

It often turns out that when we work with individual groups in 
a calculation, one or more groups ends up as a negative number. 
In the same way that we “carried” a | in melding 31129 into 
4129, we can borrow a | from a group to add 100 to the group 
to its right, as in subtraction. As an example, consider the 
calculation 54221 — 10536. We break the original numbers into 
two-digit groupings (5142121 and 1105136) and perform the 
subtraction on each piece separately, arriving at 4137!-15 = 
41361(100 — 15) = 4136185 = 43685. 

We will occasionally encounter a situation where we want to 
work with three-digit, or thousands, groupings in a number. To 
distinguish three-digit groupings, we can employ the comma as the 
separator, since commas generally separate thousands groupings 
with numbers anyway. The final melding process reduces the 
number to positive three-digit groupings by carrying or borrowing 
as before from groups to the left. For example, 217,39,2820 = 
217,039,2820 = 217,041,820 = 217041820. The group 39 was 
expanded to 039 to fill its grouping; remember, these are separate 
three-digit numbers until they are melded: 


ai 
039 
2820 
217041820 


Rarely, we will need to indicate single-digit groupings. | have 
chosen a double vertical bar (|! |) as the separator. To illustrate, 
we will show later that multiplying a two-digit number ab by 11 
gives a result al | (a + b)! |b, so 


11 x 78 = 711(7 + 8)118 = 71115118 = 8115118 = 858 
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It quickly becomes very easy to meld groups into a final answer 
in a left-to-right manner as the answer is being read out. Do you 
see that 211451-35141-3 = 344650397? 

The examples in this book have been chosen to provide realism 
to the presentation. Although it is impossible to select numbers 
that do not have some sort of special property in an application, 
I have endeavored to select values that are not particularly special 
(but, ironically, are then atypical). In a few cases | have deliber- 
ately chosen a value demonstrating the particular advantage of a 
certain technique. Where a special characteristic of a number does 
appear, as in that case, | have attempted to point it out. I have 
no interest in distorting the efficacy of the methods in this book. 
In practice we usually find that shortcuts abound. 

It may appear in some instances that the results are taken to 
extreme accuracy. However, it is my experience that results to 
lesser accuracy become easy only when ones to higher accuracy are 
attempted. Too, a good deal of personal enjoyment can be had in 
pursuing a calculation to our best efforts. Finally, this enables us 
to judge the capabilities of various techniques. 

A mathematics professor of mine once told me that he was 
surprised occasionally in calculus, rarely in trigonometry, and 
never in algebra. | wonder if he would be surprised by some of the 
arithmetic in this book. 


Bibliography 


1. James Russell Lowell, Among My Books Vol. I, Houghton, 
Mifflin and Co., Boston, 1904, p. 182. 

2. Steven B. Smith, The Great Mental Calculators, Columbia 
University Press, New York, 1983. 

3. Douglas R. Hofstadter, Gédel, Escher, Bach: An Eternal Golden 
Braid, Basic Books, New York, 1979, p. 567. 


Primitives 


The whole of Mathematics consists in the organization 
of a series of aids to the imagination in the process 
of reasoning. 


Alfred North Whitehead (1898) [1] 


<4 

Primitive” methods are those strategies for performing arith- 
metic operations such as addition, subtraction, multiplication, and 
division, as well as the means of error checking the results. These 
have been culled from a number of sources, some of which provide 
them as particular “tricks” without an appreciation of the exten- 
sions necessary for developing general tools. My intention is to 
describe methods useful in a variety of calculations; excluded are 
tricks specific to certain numbers (except for divisibility tests) or 
uncommon digit distributions within numbers. The techniques in 
this chapter, aside from their intrinsic value, form helpful tools for 
mentally calculating the functions elsewhere in the book. 


Addition and Subtraction 


Addition and subtraction offer little, actually, in shortcuts. In 
general, it is a great help to mentally group digits into twos or 
threes for the same reason that it is easier to add 49.55 to 39.62 
than 4,955 to 3,962. As mentioned in Chapter 1, the notation 


9 
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alb denotes the digits a followed by the digits b, where b is 
reduced to two digits by melding higher digits into a. In the 
above example, 


4955 + 3962 = 49155 + 39162 
= 881117 
= 89117 
= 8917 


Also, 4955 — 3962 = 101-7 = 9193 = 993. The usefulness of this 
notation will become apparent later. 
It is often convenient to adjust the numbers for simplicity: 


688 + 443 = 700 + 443 - 12 

1514 — 688 = 1514 — 700 + 12 
Some authors contend that one should never subtract per se, 
but rather think of adding a number to the subtrahend to get 


the minuend on a digit by digit, or in our case a group by 
group, basis. 


4955 Se eS 
— 131942 6 + 9 = 15, carry 1 to subtrahend 
993 9 + 1 + 9 = 19, carry 1 to subtrahend 
aoe + © = 4 


Or, 
39162 + 101-7 = 49155 


What is apparent here is that these methods all boil down to 
the same thing in slightly different ways. Believe me, though, it 
is the way we think about the numbers that can make them 
seem easy to calculate. This is particularly true in the multi- 
plication routines. 


Multiplication 


Splitting numbers into groups, or using “group vision” as 
Menninger [2] calls it, is also useful in multiplications. For 
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example, to multiply a number by 5 we halve ten times the 
number. This halving process is eased if we advantageously split 
up the number into groups of even numbers: 

5629432 x 5 = 56 2 94 320+2 
= 28147160 


Mental calculators often performed multiplication as a series of 
additions according to the Distributive Law of Algebra: 


386 x 471 = 300 © 400 + 300 © 70 + 300 #1 
+ 80 e 400 + 80 © 70 + 80 ® 1 
+ 6° 490-4 6° 70 +68 1 
with appropriate shortcuts thrown in (e.g., any calculator would 
complete two-digit multiplications almost instantly, shortening 
the steps). A digit by digit solution is possible by performing the 
classical multiplication routine vertically instead of horizontally, 
called cross multiplication [2]. To find 386 x 471: 
6°1=6 
Ge / + Srenl = Oncame5 
6°478e¢7+3¢1+5 = 8, carry 8 
8e4+3¢7+ 8 = 1, carry 6 
3°©4+6=18 
Answer: 181806 
While this eliminates the need to remember the intermediate 
values in the classical algorithm, the left-to-right nature of the 
first technique is better suited to writing the answer while it is 
being obtained. At least one contemporary lightning calculator 


uses two-digit by two-digit multiplications in cross multiplica- 
tion [3]. 


12 Dead Reckoning: Calculating Without Instruments 


For those of us without the practiced rhythm necessary for 
quickly doing these steps, there are several ways of simplifying the 
calculation by manipulating the numbers into simpler ones. In 
fact, this is where the fun lies in all of this. 

One way of simplifying a multiplication is to make use of 
the identity 


(a+ b) © (a +c) = a(a+ b+ c) + be (1) 
Therefore, 

105 x 117 = Beerer ee 5 © 17 

98 x89 = 87 e@0 + 2°¢ 11 
85 x 117 =TO2"es00 15 © 17 
= 10200 - (22 e 10 + 5 @ 7) 
or 
= 10200 —- (12 © 20 + 3 e 5) 

It is impossible to adequately convey the usefulness and perva- 
siveness of this technique in practice; if you don’t use it, do. 
Interestingly enough, this method was once taught to school- 
children in lieu of teaching the multiplication tables above 5 x 5 


[4]. To illustrate, the product 8 x 7 was found by forming 
this diagram: 


& 7 where the 2 and 3 are differences from 
Nie 10 of 8 and 7, 5 is 7 — 2 or 8 — 3, and 
aN 6 is 3 x 2, giving 56 as the answer. 

Ss 2 6 

5 


Why not? Perhaps a second-grader today should learn this tech- 
nique instead of learning multiples of 12. 

There is an analogous method in finger-counting that utilizes 
this relation as well. Again, to find the product 8 x 7, we begin 
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with two closed fists. We then raise fingers on one hand to count 
from the first number 8 to 10, and we repeat this for the number 
7 on the other hand. Then the sum of the closed fingers is the 
tens digit and the product of the raised fingers is the units digit. 

We can also extend Equation 1 to cases where both numbers 
are near numbers of which one is a multiple n of the other. A 
straightforward approach would be to multiply the smaller number 
by n, use Equation 1, and then divide the result by n. An easier 
way makes use of the identity 


(a + b) © (na + c) = al(na + c) + nb] + be 
36 x 157 = 40(157 -4 © 4) + 4 © 3 
= 40 ¢ 141 + 12 
= 5652 
Again, b and ¢ are signed quantities. 
The last example also demonstrates a convenient tool, useful 
when a multiplier is divisible by 9: 
36 x 157 = (40 — 40/10) © 157 
= (1 — 1/10) © 40 @ 157 
= 62180 —- 6128 
= 5652 


Of course, this is a special case of the common practice of 
multiplying by a nearby round number and correcting the result: 


26 x 87 = 26 © 90 - 3 @ 26 = 23140 — 78 = 2262 


When multiplying two numbers a5 and b5 that end in 5, we 
might remember the relation, 


b 
) 12s 


(a5) © (b5) = (xb > 
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35 x 75 = (Zier 
= 2625 


If (a + b) is odd, the .5 emerging from the average is assumed to 
merge into the | 25, producing | 75. 


65 x 235 = (8er 14a 
= [5775 


A final instance of algebraic manipulation to ease multiplica- 
tions is given by D. E. Smith [5]: 


(10a + a) @ (10b + b) = [(10a + a)b + ab] © 10 + ab 
22 x 44 = (22° 4+2¢4)¢10+20¢4 
= 960 + 8 = 968 
Other points to be aware of include 
25x n=nl00 + 4 
75 x n = 3n!l00 + 4 
125 x n = 8n,000 + 8 
15 & fee) © 10) or (n/2) @ 3 © 10 
11 :eeal = ane + -b) I Ib 
Here b is the units digit of a number, and a can contain one or 
more digits. Considering the last relation above, if we multiply 
al |b by a multiple of 11, we can multiply al |b by the multiple 
and continue as if multiplying by 11. I reluctantly include this 


rule, as 11m = 10n + 1 anyway, but it may be useful yet. To repeat 
our earlier example, 


22 x 44 — 11 x 88 = 81116118 = 968 
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Notice that the middle digit is melded to become a single digit, 
as indicated by the double bar notation. For larger numbers, 


44 x 643 > 11 x 25172 


= 2LIaV9 
ns Bh (Oi 2 
ee CU 


This may be extended to include multiplication by, say, 111: 
111 x al ib =al[(a +b) @ 11JIb 


where the middle section is now melded to two digits. 

Note that often we can factor a number into one or more of 
these convenient numbers. Conversely, for example, the number 
37 can be represented as 111/4. 

Squaring is an extremely important operation to perform men- 
tally, mainly because ordinary multiplications can usually be 
simplified into a square and a much simpler multiplication. For- 
tunately, squares offer advantages to us. 

Squares can be obtained by using Equation | to reduce the 
number of significant digits: 


(32)? = 30 ¢ 34 + (2) 
(185)? = 200 © 170 + (15)? = 34000 + 225 
= 34225 


Numbers ending in 5 can be squared by a trivial result of Equa- 
tion 1, often presented as a number trick: 


(35)? = [3(3 + 1)}125 
= 1225 
Conversely, general multiplication can make use of squares as 


follows (called the rule of quarter squares and of very ancient 
origins) [5]: 
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38 x 32 = (35) 

62 x 74 = (68)? — (6)? = 70 © 66 + 4 - 36 
It’s useful to me in getting the signs right to remember that the 
average number squared is larger than the individual numbers 
multiplied. 

Another convenient tool is to use the relation 

(a + b)* = a2 + bla + (a + b)] 

(31)2 = (30)2 + (30 + 31) 
(58)? = (60)? — 2(60 + 58) 

The square of a multiple of 3 can be found by extracting 
(3)? = 9 out of the quantity squared and letting this equal 
(10 — 1) aah 

(3a)? = 10a2 — a? 
(132)* = 10(44)? — (44)? 
= 1193160 — 19136 
= 17424 

One of the few explicit “tricks” I find worth remembering is 
extremely useful for squaring numbers from about 30 to 70. The 
result is obtained by adding the difference b from 50 to 25 and 
appending b’, melded to two digits [2]. 

(58)2 = (8 + 25)164 = 3264 

(37)* = 121169 = 1369 
For squaring a number between about 400 and 600, add the 


difference from 500 to 250 and append the square of the last two 
digits melded to three digits: 
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(531)? = 281,961 
(564)? = 314,391 196 
= 318,096 


Also, it is sometimes convenient to recall that when squaring 
a two-digit number ending in 1, 


(al 11)? =a7/I2all1 


where the 2a term is melded to a single digit. This follows from 
the expansion of (10a + 1). In addition, when squaring a number 
ending in 25, as in al 25 [2], 


(al25)* = 10(a? + a/2), 625 
(225)? = 525,625 
This follows from the identity 
(100a + 25)? = 1000(10a? + 5a) + 625 


By this point we should be adept at mentally multiplying two- 
digit by two-digit numbers. As you may have guessed, the later 
discussions in this book assume this ability. | recommend in spare 
moments sequencing through the two-digit numbers and squaring 
each one. This actually takes little time. The purpose is not to 
memorize the results but rather to formulate general strategies for 
arriving at them quickly. Remember here that “particular” tech- 
niques are really more general than they appear, because opera- 
tions on numbers near convenient ones can use the latter and be 
adjusted at the end for the difference. I think you will discover 
that there are many ways to skin these cats. 

To cube a number, it is generally easier not to square and 
multiply, but rather to use the expansion 


(a + b)> = a? + b? + 3ab(a + b) 
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For example, 
(23)? = 8027 + 180 © 23 
= 8027 + 10(21 © 20-2 @ 3) 
= 8027 + 4140 
= 12167 
Failing outright simplification of a large multiplication, we can 
make use of the following identity, apparently first introduced in 
a more complicated binary format in 1962 [6,7]: 
alb x cld = ac © 10% + [ac + bd - (a- b) © (c — a] ae 
+ bd 
6823 x 4519 = 68 © 45 x 104 + (68 © 45 + 23 @ 19 
— 45 © 26) x 10° +73 «1 


=  §060100))00 
+. 73 27100 
+ 4|37 


3083 131137 = 30833137 


thereby getting an eight-digit result by evaluating the products of 
only three two-digit by two-digit numbers. The usual expansion 
requires four such products: 


alb x cld = ac © 10* + (ad + bc) © 10% + bd 


Both equations become identical when squaring numbers. 

Knuth [6] comments that this method, while apparently never 
used by lightning calculators, should make eight-digit by eight- 
digit mental calculations “reasonable” by nesting the procedure. 
Although it reduces the number of two-digit multiplications from 
16 to 9, it does require temporary storage of a few numbers simul- 
taneously, a trivial task for a machine but a difficulty for mortals. 
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Division reveals difficulties in simplification, serving to deepen 
my suspicion that addition, subtraction, and multiplication are 
artifacts of God, while division is a construct of man. Nonetheless, 
some things can be done. Following an overview of the useful 
properties of repeating decimals, we will delve into strategies for 
reciprocals of one-digit and two-digit numbers and divisors near a 
multiple of 10 before tackling large divisors. 

To begin, we discuss certain properties of repeating (sometimes 
called recurring) decimals, that is, those numbers with digits 
(usually lying to the right of the decimal point) that recur in repeated 
groups [8-10]. For example, 1/7 = 0.14285714285714..., 
demonstrating this property. It happens that for a fraction s/t, a 
repeating decimal group will occur if t has one or more prime 
factors other than 2 or 5. Otherwise, since 2 and 5 are prime 
factors of 10, an exact decimal value will result (e.g., 13/16 = 
0.8125). A prime number is a positive integer other than | that 
is divisible only by itself and 1, such as 2, 3, 5, 7, 11, 13, ete. 
When repeating groups do occur, they may follow initial non- 
repeating digits. However, awareness of the existence of these 
groups can save a good deal of effort. 

A pertinent question at this point is how we determine the 
number of decimal places we need to either terminate the number 
completely or complete the first repeating group and, in the latter 
case, the number of digits in the repeating group. For a denomi- 
nator t with prime factors of 2 and 5 only, the number of decimal 
places will equal the highest power of 2 or 5 contained in t. In 
our example above, 16 = 24 © 5°, so the ratio 13/16 terminates 
after the fourth decimal place. Some would say that a terminating 
decimal is actually a repeating one, spewing out zeros or nines to 
infinity. We will not label these as such, if only to simplify 
our terminology. 

Unfortunately, there is no general rule for determining the 
recurring period of repeating decimals, although there are some 
useful properties of some cases. First, we can see that a group will 
begin repeating when, in long division, we arrive at a remainder 
that is the same as earlier in the division, whence the whole 
process starts repeating. Since there are only (t — 1) different 
remainders possible for a divisor t, the repeating group cannot 
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have a period greater than (t — 1). Our first example of t = 7 
illustrates a repeating group of length (t — 1), as do such numbers 
as 17, 19, 23, 29, 47, 59, 61, and 97. The reciprocal of 97 is 
occasionally asked of lightning calculators [11]. 

It turns out that for reciprocals 1/t with t prime, the period of 
the repeating group, if not equal to (t — 1), must be exactly 
divisible into (t — 1). For example, for the fraction 1/43, the group 
period must equal 42, 21, 14, 7, 6, 3, or 2. For a prime t between 
3 and 487, the period of 1/(t") equals t™! times that for 1/t. 
Therefore, the repeating group of 1/49 contains 7 © 6 or 42 digits. 
The period of a fraction s/(3") equals (3"~*) digits. 

For a reciprocal 1/t with t composite (nonprime, or factorable), 
we can sometimes break t into prime factors 22 ¢ 55 © t, © t) e 
t; ©..., where tj, ty, t;,... are all different and are greater than 
5. Then the period of the repeating group will equal the lowest 
common multiple of the periods of the repeating groups for 1/t,, 
1/t), 1/t3, etc. Yielding to a heightened sense of apprehension, | 
offer the example of 1/73920. Here we may find by factorization 
methods discussed later that t = 2© ¢ 5 © 3 © 7 © 11. The recurring 
periods of the last three factors are 1, 6, and 2, whose lowest 
common multiple 6 is the recurring period of 1/73920 = 
0.000013528138528138528 ..., where we notice that the recur- 
ring group appears after six nonrepeating digits, 000013. 

In addition, the recurring period for numbers which are factors 
of 9, 99, 999, etc. will equal the number of nines in the lowest 
such multiple. For example, 999,999 = 3 e3e3e¢7e11 e 13 
© 37 yields 63 combinations as factors. Except for the six that are 
factors of 9, 99 or 999 (3, 11, 33, 37, 111, and 333), all have a 
recurring period of six digits, including the factors 7, 13, 21, 39, 
63, 77, 91, and so on. Since 3 is a factor of 9, its reciprocal has 
a recurring period of one digit. The numbers 11 and 33 produce 
two-digit periods, and 37, 111, and 333 yield three-digit periods. 

In another special case of a reciprocal 1/t where t = 23 @ 55 e ¢? 
and t; is any prime between 5 and 487, the period of the repeating 
group equals that of 1/t; multiplied by t?-!. For 1/968, we have t 
= 23 e 11? with a recurring period of 2 © 11 = 22 digits. 

Now, for irreducible fractions s/t yielding groups of even periods 
and with t prime, corresponding digits of the first and second 
halves of the period add up to 9. Therefore, if we know that 1/7 
yielded an even period (as it does since it is one of the ones noted 
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for periods of length (t — 1) earlier), we need only find the first 
three digits .142 to know the last three of the group as 857. Again, 
3/7 = .428571.... Of course, this halves the effort required to 
find or memorize the reciprocal of, say, 97 as well. Other fractions 
with even periods but with t nonprime often exhibit this phenom- 
enon, but not always. They always will if no factor of t divides 
(10° — 1), for n a positive integer. If t is not divisible by 2, 3, or 
5 and 1/t has a period of (t — 1) digits, this situation will occur 
as well. Finally, for t = a? © b? © ..., where each multiplier is a 
distinct prime other than 2 or 5, halves will be complementary if 
and only if the periods for 1/a, 1/b, etc. contain the same power 
of 2. 

Also, for any reciprocal 1/t, the last digit of the repeating group 
will be 1 or 9 if the units digit of t is 9 or 1, respectively. If t is 
prime and the units digit is not 1 or 9, the last digit of the 
repeating group will be the same as the units digit of t. Therefore, 
the repeating group for 1/7 must end in 7, so if we take the 
fraction to two digits, or .14, we immediately know that 1/7 = 
MezOD/.... 

We turn now to the task of finding the decimal expansion of 
a reciprocal 1/t. The process of long division, which we can end 
when a remainder is repeated in the course of the calculation, is 
the obvious approach. There is another method which is fre- 
quently more convenient, particularly if the period of the recur- 
ring group is large. We divide by long division until some digits 
are determined and a low remainder is found. Then we can 
multiply the entire quotient (including the fraction remainder/t) 
by the remainder, converting the fractional part to proper form. 
Any digits to the left of the decimal point are discarded. This 
process of multiplication can be repeated until the recurring group 
is found, with checks provided by our earlier observations. For 
example, let us consider the reciprocal of 43: 


1/43 = .023'Y, 
11/43 = 11(.023'%,) = .253 121/43 = .255 35/43 
117/43 = 11(.255°%,) = 2.805 385/43 = 2.813 41/43 


C&C. 
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Then 1/43 = 0.023255813.... 

We haven’t reached the end of a recurring group yet, and of 
course it may take 42 digits. We could have continued our initial 
long division further and found 1/43 = 0.0232558 6/43. On the 
other hand, we can continue to 1/43 = 0.023255813 41/43, either 
by long division or even better as the intermediate result from 
multiplying by elevens in the calculation given above. Here we 
can take our remainder as —2 (I find that authors of arithmetic 
methods often don’t take arithmetic literally enough). Then 


1/43 = 023255814 -2/43 
~2/43 = -2(.023255814 -2/43) = -.046511628 +4/43 
~27/43 = -2(.046511628%,;) = +.093023256 -8/43 


See the recurrence beginning in the last result, realizing that 


6 — %, = 5 °% 3? We find then that 
1/43 = .023255813953488372093 


where the raised line indicates the recurring group. 

If we don’t care to work with the fractional parts at all, we can 
ignore them if we keep the digits to the left of the decimal point 
in each step [9]. These extra digits are then melded into the digits 
previously found. Again, 1/43 = 0.023 with a remainder of 11, 
which multiplies the numbers in each line: 


023 
D5 
2805 
023759005 ... 


Here, however, we would not discover the more convenient 
remainder of —2 which, as we have seen, would greatly simplify 
finding the rest of the repeating group. For reciprocals of numbers 
such as 37, 27, 7, 11, or 13, for example, that have multiples one 
less or more than a power of 10, the remainder after division to 
this number of places will obviously be 1 or —1, respectively. Then 


since 37 © 27 = 999 and 7 e 11 © 13 = 1001, we know that 
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1/37 = .027 
1/27 = .037 


and since 11 @ 13 = 143, 7 © 13 = 91, and 7 © 11 = 77, 


we find 


143 091 
= 143 = 091 
+ 143 + 091 
etc. eke: 
1/7 = .142857 1/11 = .09 
077 
_ 077 
+ 077 
GEC. 


1/13 = .076913 


Now 43 is prime and we have found above that the period of 
the recurring group is 21, a factor of (43 — 1). For t’s containing 
more digits, we can, with some luck in remainders, save a great 
deal of time with these methods. 

For numbers whose reciprocals yield recurring groups of length 
(t — 1), values of s/t with s > 1 contain the same recurring group 
with its digits cycled to perhaps a different starting position. This 
happens to be true for 1/(7)? as well. 

For example, we recall that 1/17 has a recurrent period of 
(t — 1) and we can find its value through the previous techniques. 
We need only find the first eight digits, as the last eight 
are complementary. 


1/17 = .0588 4/17 

4/17 = 4(.0588%,) = .2352 16/17 
Then 

1/17 = 0.0588235294117647 
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Here we also see that since t is prime with a units digit of 7, the 
last digit of the recurring group is also 7. Realizing this, we actually 
only needed to ascertain the first seven digits of the reciprocal. 

Now, say we’re interested in the value of 13/17. We found the 
decimal expansion of 1/17 fairly easily, so let’s divide 13 by 17 to 
a couple of digits or multiply the first few digits of 1/17 (which 
we round off) by 13, giving 13 © 059 = 767. This implies that 
the recurring group begins with the digits 764 in 1/17, and we 
assert that 


13/17 = .764705882352941 1 


which is indeed true. Lightning calculators of the showmanship 
variety occasionally use a trick involving a reciprocal such as 1/7 
to perform seemingly unbelievable multiplications. Through a 
confederate in the audience or some other contrived means, the 
first number is produced, in this case, 142857. Now, we know that 
1/7 = .149857, so 142857 = 10%1/7) - 1/77 = 1c =e 
Therefore, any other number easily multiplies this artifical one. 
For example, a legitimate audience member volunteers the 
number 6297 as the multiplier. Then, 


6297000000 — 6297 
i 


6797 « 142657 = 


_ 6296993703 
7 


and the short division is also done mentally to arrive at the correct 
answer of 899570529. Beware. 

Returning to our discussion, let us consider the general case of 
finding the decimal expansion of s/t, where s is not necessarily 1. 
First, we can round off the divisor in many instances by adjusting 
the remainder or dividend in each step of the usual division 
process. For evaluating 1241/78 = 15.910... , for example, we can 
replace the divisor by 80 and correct each remainder by the 
errant amount: 
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124/80 =1 ~~ remainder = 44 + 1 © 2 = 46 
461/80 = 5 remainder = 61 + 5 @ 2 = 71 
710/80 = 9 remainder = -10+9e2=8 
80/80 = 1 remainder =0+1¢°2=2 
etc. 


Notice that in the third step we had to back up and increase 
the previous quotient digit 8 to 9 because the adjusted remainder 
ended up greater than the divisor 80. The frequency of this 
situation increases as the original and rounded numbers are more 
distant. For an original number ending in 9, this situation never 
occurs during the string of zeros following the significant digits in 
the dividend. Digits in the quotient are simply inserted (added to 
zero) into the dividend a number of places to the right given by 
the power of 10. In this case, it is convenient to divide both the 
divisor ending in n nines (adjusted to a number ending in n zeros) 
and the dividend by 10n. Then as the division progresses, each 
digit of the quotient is added to the dividend n places later. For 
example, in the reciprocal 1/29, we round 29 to 30 and divide this 
and the dividend by 10: 


03448275... 
3).103448275 ... 


where we add each digit of the quotient, as it’s obtained, to the 
dividend one place later. 

Now in the situation where a fraction s/t has a denominator t 
of the form (10m — 1) for m an integer, that is, it ends in a 9, 
we can mechanize the previous method of division to speed the 
process [9]. We can think in terms of the following relations: 


ag = s (not necessarily 1) 

b,, = integer part of a,/m 

C, = remainder of above 
an+t = 10c,, ” by 
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Now b, will consist of the digits of the quotient. The term c,, 
the remainder when a, is divided by m, may be written as 
a, mod m, where c, is said to be congruent to a, with respect to 
the modulus m. Let us consider the fraction 1/29 again. Then 
m = 3 and we have the following sequence of steps: 


— 

Nl 0 | 
10 3 | 
13 4 1 
14 + 4 
24 8 0 

8 2 2 
22 7 1 
17 5 ] 
etc. 


giving us as we ripple mentally through the rhythmic process: 
1/29 = .03448275... 
For the fraction 19/29 we have 


a ae 
19 6 
16 5 
15 5 
5 1 
21 7 


oNnore|o 


Gtc. 


correctly giving 19/29 = 0.65517.... 

Any divisor ending in 3 will transform into one ending in 9 
through multiplication of the numerator and denominator by 3. 

For t of the form (10m + 1) or in other words ending in 1, and 
for s not necessarily unity, we set a,,; = 10c,, — b,. However, since 
we can end up with negative values of b, (which may be melded 
to obtain a correct result), we can ensure positive values by setting 
b,, to one less than a,/m if m divides a, without remainder. For 


1/21 we have 
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a  b & 
1 0 1 
10 4 2 
16 7 2 
13 6 1 
4 1 2 
19 9 1 
1 0 1 


and here we have the initial line repeated, revealing that 
1/21 = .047619 


Notice that divisors ending in 7 can be converted to ones ending 
in 1 by multiplying the numerator and denominator by 3. 

If we don’t mind obtaining the digits in a right-to-left manner, 
we can use another method, given by Taylor [9]. For a denomi- 
nator of the form (10m — 1), we can add shifted terms multiplied 
by m, beginning with the numerator itself. Taylor gives the 
examples of 


1/39 = 025641 and 1/13 = 3/39 = .076923 
1 
4 3 
16 12 
64 48 
256 192 
1024 768 
+ 4096 + 3042 
. . 025641 ... OMe 


Taylor gives other methods suitable for very particular divisors 
(such as squares) that we will not consider here. 

For divisors generally lying near a power of 10, the remainder 
can be also be adjusted while the entire dividend is divided by the 
power of 10 [2]. For example, for 129641/97 we can perform the 
division by 100 first: 


28 Dead Reckoning: Calculating Without Instruments 


129641/100 = 1296 remainder = 41 
new remainder = 3 @ 1296 + 41 = 3929 


3929/100 = 39 remainder = 29 
new remainder = 3 @ 39 + 29 = 146 
146/100 = 1 remainder = 46 


new remainder = 3 @ 1 + 46 = 49 


Therefore, 129641/97 = (1296 + 39 + 1) = 1336 with a remainder 
of 49. Of course the division process can be carried further into 
the region to the right of the decimal point. 

For divisors that have a multiple that lies near a power of 10, 
the remainder can also be adjusted as in the previous method. The 
difference consists of multiplying the quotient in each step by the 
value of the multiple factor. To illustrate this, we can find the 
value of 4330463/332, where 1000 = 3 © 332 + 4. 


4330463/1000 = 4330 remainder = 463 
new quotient = 3 © 4330 = 12990 
new remainder = 4 @ 4330 + 463 = 17783 
17783/1000 = 17 remainder = 783 
new quotient = 3 @ 17 = 51 
new remainder = 4 @ 17 + 783 = 851 
851/332 = 2 remainder = 187 


Therefore, 4330463/332 = (12990 + 51 + 2) = 13043 with a 
remainder of 187. 

Now let us consider a method of cross division, or Fourier 
Division, as proposed by Joseph Fourier (1768-1830), useful for 
calculations involving multidigit divisors [12-14]. If we are divid- 


ing C by A to get B, then AB = C and 


A= a,la,lazlay... 
x= b, | bz 1b; |! by... 
Pil qi 
P2!qQ2 
p3!q3 
p4! qa 


©: = -e; bepl'agliae.. . . 
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where the partial products p,|q, are formed by multiplications of 
a and b terms. In particular, 


Pl lq = a,b; 
p2!q2 = ayb2 + azb, 
p3!q3 = ayb3 + ayb7 + ab, 


or, 
n 
Pr! dn = y a by+1-k ms ©, 2 3,... 
k=1 


where the Greek letter } (sigma) indicates the sum of the terms 
for k from 1 to n here. Also, 


lc, = Ip, + |q,_1 with qo = 0 


We are interested in solving for B = b,!b,|b;... given A 
and C: 


_ Pil a _ alle, — p») 


ay ay 


by 


‘eile, 


with a remainder R, = | p2 
ay 


-. Pa! qa —eagb; Ryle; — p;) —-azb, 
eS eee ae 
a| ay 

a R, 1 c3 = arb, 


;R2= Ip; 
aj 
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oe P3!q3 — arb) — asb, 


b; 
aj 
_ Ral cg — p4) — agb2 — aby 
a) 
Ry I c4 = arb = a3b, 
= ER; = Ip, 
ay 
In general, 
n 
Pa | a5 = bE Anby+1-k 
E—2 
b, = 
a} 
n 


Pa! (casi = Diwat a a ae 


aj 
nN 
Pn! (essai 9 = > Anbn+1_k 
F=2 Pn+1 
a} a} 


and, defining the remainder R,, as | p,.1, 


ig 
Replicas — y Anba+1-k 
k=2 
5, = — °° °° = = I Pati 
1 


As an example, consider the division of C = 42472482 by 
A = 874921. We consider two-digit operations here, although 
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single-digit operations are possible after, say, a two-digit a, is split 
off as the divisor. Notice that for accuracy the remainders R,, are 
adjusted to give numerators which are less than a,/2, occasionally 
resulting in negative values of b. 


C 42147124182 
A 87149121 


42147 
b, = —— = 49; R, = -16 
87 


quotient thus far: 49 


_ 16124 - 49 © 49-1576 - 49 © 49 


b2 = 
87 87 
= — (46,R -25) = — 46; R; = 25 
49 |-46 
25182 — (-46) © 49 - 49 @ 21 
ae  C_ 44. R, = 21 
87 
49 |_-46 | 44 
21100 — 44 @ 49 + 46 e 2] 
ee 6 2) oe 165 
87 
= -38 ; Rg = 16 
49 |-461 441-38 
16100 + 38 @ 49 — 7 | 
ee e499. R, = 15 


87 
49 |-46 1441-38129 
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Melding the b’s and fixing the decimal point results in B = 

48.54436229.... We can continue the process to additional 
decimal places if desired. Notice that there is no increase in 
difficulty for steps after b), since A has only six digits. Nonethe- 
less, a certain knack has to be acquired to perform the algorithm 
with any speed. 
If the first two digits of C are larger than the first two digits of 
A, then c,!c2/a; will produce a three-digit (and therefore inac- 
curate) quotient. In these cases, we can precede the first digit of 
C with a zero (e.g., C = 92472482 becomes 091241721481 20). 
In addition, we should round a, to the nearest integer; if A were 
to equal 876221 we would transform A into 88!-—38121. Finally, 
it is again worth mentioning that sometimes division by a partic- 
ular number can be eased, as in multiplication. For example, 
a/l5 = 2af(3 @ TQ). ,etc. 

Since we occasionally encounter the situation where we precede 
C with a zero, perhaps we should add another example where a 
small number is divided by a four-digit number, say 23/1024 = 
0.022460938 .... We will add a zero to the end of C as well to 
fill out the two-digit quantity left hanging when the zero is added 
to the front; we'll sort out the decimal place at the end. 


C = 02130 
A = 10124 
B = C/A = b, |b) !b3... 


02130 
b,; = —— = 23; R, = 0 
10 


23 
00100 - 23 © 24 
a = — (55)R 2) = 55) Rp = 
10 
231 -55 
-2100 — (-55) © 24 
a .- 112; =O 


10 
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Since b; is three digits long, let b) = —54 and R, = —12. 
Then 


-12100 - (-54) © 24 
eee ate 8.9. R, = 4 


10 
23154110 
-4100 — 10 © 24 
i ee 6, By = 0 
10 
231-541101-64 


Melding the b’s and adjusting the decimal point yields, without 
too large a total effort, the result B = 0.022460936..., and the 
next step would correct the inaccuracy in the last digit. 

While cross division and other methods can be used in calculat- 
ing a reciprocal exactly, other methods may prove convenient for 
finding an approximation to the reciprocal of a number near a 
round number. 

For approximate reciprocals of multidigit numbers, we can use 
an iterative relation for 1/t [15]: 


as =e (2 - tx,) (2) 


This relation is derived by assuming that an approximate value 
x, is at hand. Then the quantity (x, — 1/t) is minimized by 
minimizing (1 — tx,,). Now to ensure that (1 — tx,,;) < (1 — tx,), 
we write 


1 — tx,4, = k(1 — tx,)? 


For k = 1, we arrive at Equation 2 by solving this equation for x,,,1. 
The process is a second-order one, doubling the number of 
correct digits in each iteration. For t = (a — b), xo = I/a, and with 
(signed) b small relative to the round number a, we can write the 
formula for x; in a form more suitable for mental calculation: 


X} = Xo + ba? 
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where we can divide b by a twice rather than by a’ once to 
simplify the division. We can take as an example the reciprocal 
of 387, or 0.00258398 ...: 


a = 400 
b = 13 
xo = .0025 


0025 + 13/(400)? 
0025 + .0325/400 
00258125 


xX] 


We would generally not calculate x, because of the number of 
digits in the numbers involved. However, we can get greater 
accuracy if we choose to use a third-order relation [6], 


Xena) = Kyo Re ee xe! — tx.) 
or, 

Keo, = Xalle+ (1 Steal + (1 —D] 
or, in the earlier notation, 


b + b2/a 


a 


X; =e + 


Repeating the example above, 


13 + 169/400 


x, = .0025 + 
400 © 400 
13.4225 
= (05 + = = 
400 © 400 


00258389 ... compared to .00258398 . .. 
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Finally, it may be useful at times to convert from a repeating 
decimal to a fraction. We can limit ourselves to one method as 
follows. We can see that for a reciprocal 1/t yielding a recurring 
group g of period p, t must divide into (10° — 1), since it is this 
remainder of | that starts the whole process again. In fact, for the 
fraction s/t, 


= = 
10° ~ 1 


S 
Ct 


In the chapter on extracting roots, we will see that the lightning 
calculator A. C. Aitken at least implicitly used this technique. 
From an initial guess of 50/7 = 7.142286... for the square root 
of 51, he found a closer value of 7.141429. Therefore, he took as 
a new initial estimate the value 7.141414... , which he converted 
into a convenient fraction: 


-_ 14 
om = 7 + ———— =707/99 
107 - 1 


The Greatest Common Divisor 


It is also occasionally useful to find the greatest common divisor 
of two numbers a and b, or gcd(a,b). This is especially useful for 
simplifying the division of a large number by another large 
number, and it is also utilized in other calculational techniques. 

The traditional method of finding the greatest common divisor 
of two numbers, Euclid’s Algorithm, dates from 300 — 400 B.C. 
[6,7]. The process is iterative and works for any positive integers 
a and b. In each step, the larger number is replaced by the 
remainder when the larger number is divided by the smaller. 

As mentioned earlier, we write the relationship between num- 
bers having the same remainder when divided by b as c =a mod b 
for a > b; in the terminology of congruences, ¢ is congruent to a 
modulo b. When finding the remainder r when a is divided by b, 
we will use the traditional equals sign (a mod b = r). The 
reasoning behind Euclid’s Algorithm lies in the rules of modular 
arithmetic, which allow straightforward addition, subtraction and 
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multiplication of terms of this form. For example, since 230/4 
leaves a remainder of 2, we write the congruence 2 = 230 mod 4. 
By simple modular arithmetic, we can also write that 2 = 200 mod 4 
+ 30 mod 4, 2 = 292 mod 4 — 62 mod 4, and so forth. Modular 
multiplication is just as simple; the congruence 2 = 230 mod 4 can 
also be expressed as 2 = (10 mod 4) x (23 mod 4). In terms of 
remainders, we have 230 mod 4 = 2, (200 mod 4 + 30 mod 4) 
mod 4 = 2, (292 mod 4 — 62 mod 4) mod 4 = 2, and [(10 mod 4) 
x (23 mod 4)] mod 4 = 2. These are easily verified. Remember that 
if we end up with a negative result, we add 4 until we arrive at 
a positive number. From above, (292 mod 4 — 62 mod 4) mod 4 
= (0 — 2) mod 4 = —2 mod 4 = 2 mod 4. 


Division is more complicated, for if cd = ad mod b, then 


c =a mod 


gcd(d,b) 


At any rate, to exercise Euclid’s Algorithm let us simplify the 
fraction 2745/13664: 


13664 = 4 © 2745 + 2684 
2745 = 1 © 2684 + 61 
2684 = 44 © 61 +0 


Therefore, gcd(13664,2745) = 61 and we can transform the 
division to 45/224. 

We can modify this procedure by choosing the multiplier of b 
in each step to obtain the minimum absolute value of the 
remainder. This sometimes gives us a negative remainder, of 
which the absolute value is used in the next step. 


13664 = 5 ¢ 2745 - 61 
2745 = 45 ¢ 61 + 0 


This is called the least-remainder algorithm for finding the great- 
est common divisor of two integers [16], and the number of steps 
is decreased over Euclid’s Algorithm by the number of negative 
remainders that occur. 
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There is another way of reducing the number of steps and 
division effort generally required by these algorithms [17]. Notice 
that we might as well have divided out all multiples of 2 from the 
even value 13664, since it obviously does not divide 2745. If we 
do this and start with two odd numbers, we can divide out of the 
remainder a mod b in each step all powers of 2 until it is odd. If 
the remainder is odd to begin with, we can increment the 
multiplier of b to arrive at an even (and now negative) remainder. 
We find, for our particular example above, that there are savings 
only in the magnitude of the original numbers since the steps were 
few due to the fortunate instance of a remainder a mod b very near 
b. To illustrate, 


gcd(13664,2745) = gcd(427,2745) 
2745 =7° 427-48 61 
477 =7¢« 61+0 


giving the result gcd(13664,2745) = 61. As another example, 
compare the two methods for finding gcd(28567,3829) = 7: 


Euclid’s Algorithm Modified Euclid’s Algorithm 


28567 = 7 © 3829 + 1764 28567 = 7 © 3829 + 4 © 441 
9629 = 7 0 1164+ 301° 38629 =9 © 441-4 35 
1764 = 5 © 301 + 259 441=11¢35+8e/7 


301 = 1 © 259 + 42 er § © 7+ 0 
259 =6 ¢ 42 +7 
42=6°e7+0 


In general, there is a real savings in doing this. In addition, we 
can eliminate multiples of other primes from the remainders if 
these primes are eliminated from the original values a and b. 
Again, our first example collapses prior to demonstrating the 
effect, since extracting, say, a multiplier of the prime number 5 
from 2745 leaves 549 and 549 = 1 © 427 + 2 © 61. However, in 
the second example we can immediately see that neither 28567 
nor 3829 contains a multiple of 5. Therefore, at the second step 
in the modified Euclid’s Algorithm above, we can reduce the 
already reduced remainder 35 to 7. The operations become 
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28567 = 7 © 3829 + 4 © 44] 
3829 =9 © 441-4507 
441 = 63 ©7+0 


where we now eliminate both the 4 and the 5 in the second line. 

If we don’t care to perform the division at all to find a mod b, 
we can use a method of halving operations throughout. This 
approach was proposed in 1962 for computer applications [6]. The 
process involves initially halving a or b until they are both odd 
(at least one is assumed odd to begin with or the common factor 
of 2 is obvious). Then in each step the absolute value of the 
difference (a — b) is halved until an odd value results, which then 
replaces the larger of a and b. 

To return to our earlier example of 2745 and 13664, the even 
value 13664 is halved until a value of 427 is obtained. Then 


2745 — 427 = 2318 — 1159 
1159 — 427 = 732 — 183 
427 — 183 = 244 — 61 
183 — 61 = 122 + 61 
61-61 =0 


This algorithm is actually a manifestation of the previous 
method, with the multiplier in each step being unity and only 
multiples of 2 being dropped. It is apparent that computers can 
perform halving operations as simple single machine cycle shifts 
of binary numbers, so this is a major improvement for them. On 
the other hand, we are capable of perceiving when higher divi- 
sions can be performed, since a number is divisible by 4 if its last 
two digits are divisible by 4, and by 8 if its last three digits are 
divisible by 8 (since 4 and 8 divide 100 and 1000, respectively). 
This significantly decreases our steps in both of the above algo- 
rithms. In addition, we can notice that in the fifth step above, 183 
is a multiple of 61, so we can terminate the procedure early. 
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The lowest common multiple (lcm) of two or more numbers can 
be found from their greatest common divisor. Earlier, we determined 
that gcd(13664,2745) = 61. The number 61 then contains all 
prime factors common to both 13664 and 2745. Then 13664 = 61 
e 224 and 2745 = 61 © 45, giving lcm(13664,2745) = 61 ¢ 45 
e 224. This may be useful for adding and subtracting fractions: 


11 7 mengs + 7° 224 
i 
13664 2745 61 © 45 © 224 


_ 2063 
614880 


Error Checking and Divisibility Tests 


Error checking of calculations is well worth the little extra 
effort. A result, of course, should be automatically checked for 
reasonableness in the overall sense of correct number of digits, 
reasonable initial digits, and exact last digit. 

Other methods generally consist of divisibility tests, and are 
usually used in addition, subtraction, and multiplication opera- 
tions. Some divisibility tests provide a remainder if divisibility 
does not hold, and may be used directly as a division process for 
one- or two-digit divisors. 

Casting out nines, for example, represents a divisibility test for 
9, and casting out 99’s can improve the accuracy from 89% to 
99%. Casting out nines (the nines test), if you recall, amounts to 
finding the sum of the digits of the number and iterating this 
process to arrive at a single digit representing the remainder when 
the number is divided by 9. Intermediate digits like 9 or digits that 
sum to 9 are easily eliminated during this process. For example, 


236439 mod 9 = (2 +3 +6 + 4+ 3 +9) mod 9 
= 27 mod 9 
= (0 
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We can also eliminate the “36” combination, since 3 and 6 sum 
to 9, as well as the lone 9 at the end. The remaining digits (2, 
4, and 3) sum to 9, so our answer is 0. We can test the calculation 


236439 x 15 = 3546585 as follows: 
236439 mod 9 x 15 mod 9 = 3546585 mod 9 
ox6=0 check 


Casting out 99’s works on digits pairwise from the right. As an 
example of casting out 99’s, consider the calculation 229 x 721 


= 165109: 
2129 x TAZE =aGil SOD 
test: 31 © 28 = 76 
8168 = 76 
76 = 76 check 


The main drawback of the nines test, and one that in practice 
lowers its accuracy below 89%, is its failure to detect errors in 
place. That is, if a partial result is shifted an incorrect amount 
(894 instead of 8940, for example) before being added to the sum 
being accumulated, it will not be detected in the final result. 

A better tack is to cast out elevens, which amounts to subtract- 
ing the sum of the even-place digits from the sum of the odd-place 
digits. If a negative result appears, multiples of 11 are added to 
bring the overall result to a non-negative integer less than 11. 


965109 mod 1D=(9 + 1 +6)-(0+5+1)=10 


This elevens test detects errors of place and provides a true 91% 
accuracy. We can cast out 101’s by subtracting the sum of even 
pairs of digits from the sum of odd pairs of digits, and so forth. 
In general, one can construct divisibility tests for various 
numbers, and these are helpful not only in many division prob- 
lems, but also for factoring a number whose logarithm is sought. 
Other than the nines test and elevens test we’ve mentioned, we 
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have already noted in the last section that a number is divisible 
by 4 if its last two digits are divisible by 4, and by 8 if its last three 
digits are divisible by 8, since 4 divides 100 and 8 divides 1000. 
In addition, we see that a number is divisible by 3 if the result of 
the nines test is divisible by 3. 

One convenient method of testing divisibility by a divisor d of 
a number N is to recall a multiple m of d that lies near a power 
t of ten, the remainder designated here by r [2]. 


md = 10°+r or 10° mod md = -r 
If we dissect N into t-digit groupings n,, 

N = np + n,10' + n)10% + n310* +... 

= no + 10%(n,; + 10"(n, + 10(n; +...))) 

and therefore, 

N mod md = [ng — r(n, — r(n) — r(n3 —... )))] mod md 

Simply put, then, one works from the left end of N, subtracting 
each rn, from n,_;. For the resulting number R (consisting of at 
most t digits): 

R mod d = N mod d 
which directly provides us with the remainder if desired. This 
technique explains why the elevens test works as it does. 

For example, to cast out 17’s, we know that 6 @ 17 = 100 + 2. 
Therefore, r = 2 and since 100 = 10’, we separate N into pairs 
of digits. Beginning on the left, we subtract twice the digit pair 


from the next pair, dropping multiples of 17 from any pair when 
convenient. For N = 165109, 


16151109 + 19109 4 2109 5 5 


and 165109 leaves a remainder 5 when divided by 17. As an aside, 
we've also found here that 165109 mod 6 = 5 as well. 
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There are a number of convenient divisors. For example, divi- 
sion by the prime number 37 is eased because 37 © 27 = 1000 - 1. 
In addition, since each group can be reduced by any multiple of 
the divisor, and 3 ¢ 37 = 111, we can easily subtract multiples 
of 111 from each triplet: 


784,165,109 — 7,54,-2 — 61,-2 — 59 > 22 


Another one of interest is 7 © 11 ¢ 13 = 1001, useful for simplify- 
ing the tests for any or all of these factors: 


1,109,185 — 108,185 — 77 


Therefore, 1109185 is divisible by 7 and 11 and leaves a remainder 
of 12 when divided by 13. 

We can use a related algorithm [18-20] for determining the 
divisibility of a number N by an odd integer d whose last digit is 
not a 5. A disadvantage of this method is that, except for d = 3 
or 9, the remainder N mod d does not emerge if nonzero. 
However, the method is interesting and can be useful. 

The basic approach for a number N consisting of a sequence of 
digits is to add a multiple m of the terminal digit to the number 
formed by the remaining digits. This process is iterated, and each 
resulting number N,, N2, N3, etc. is divisible by d if the original 
number N is. We terminate the process when we easily see 
whether or not a value of N is indeed divisible by d. 

We consider a divisor d consisting of digits D,d5, where d) is 
a single digit and D, can be multidigit or single-digit, including 
O. The multiplier m is given for d, = 1 as (—D,) and for d, = 9 
as (D; + 1). The absolute value of m here represents the multiple 
of 10 that N is either one greater than or one less than; the 
negative sign for d, = 1 corresponds to the alternating signs in the 
elevens test. For d) = 3 or 7, we multiply d by 3 and arrive at one 
of the cases above. 

For example, the number N = 304 can be tested for divisibility 
by d = 19 as follows: 


N = 304 
N, = 30+2¢ 4 = 38 
N, =3+2¢8=19 
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Therefore, 304 is divisible by 19. For d, = 9 without multi- 
plication by 3, the process will always terminate in d if d divides 
N. In fact, for the number ¢ constructed of final digits of all but 
this last N, the quotient q is given by (10° — c), where t is one 
less than the number of digits in N. Here, q = 100 — 84 = 16, so 
304 = 19 © 16. 

For d, = 1 without multiplication by 3, the process will always 
terminate in zero if d divides N, and then gq = c. As an example, 


N = 83049 d=3l 
N, = 8304 -3 ¢ 9 = 8277 
N, = 827 —-3 © 7 = 806 
N; = 80-3 © 6 = 62 
N,=6-3¢2=0 
and 
q = 2679 
When d, becomes 9 or 1 upon multiplying d by 3, the process 


will end with a multiple of d if d divides N. Extraction of the 
quotient is not straightforward. As examples, consider 


d = 13, becoming 39 d = 27, becoming 81 
N = 15587 N = 2619 
N,; = 1558 + 4 © 7 = 1586 N, = 261 — 8 © 9 = 189 
N) = 158 + 4 © 6 = 182 N, = 18-8 @ 9 = -54, 
N; = 18 + 4 @ 2 = 26, divisible by 27 


divisible by 13 


An extension of this strategy consists of adding or subtracting 
a multiple of a group of terminal digits of N. While in most cases 
this is too cumbersome to be of much practical use for mental 
work, a few cases may be useful [20,21]. For example, divisibiliry 
of N by d = 29 or 23 can be determined by subtracting twice the 
number formed by the last three digits of N from the number 
formed by the remaining digits, and then iterating. For N = 5851417, 
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N, = 5851 —2 © 417 = 5017 
N,=5-2¢17=-29 


which shows divisibility by 29, but not 23. The quotient is hidden. 
The source of this technique, as you probably have guessed by 
now, lies in the fact that 2001 is divisible by 23 and 29. Based on 
this principle, Smith [3] lists the methods lightning calculator 
Wim Klein might use to determine possible factors of the number 
N = 114043 (note that these methods don’t directly provide 
correct remainders if nonzero). In each case shown in Table 1 the 
reduced number (and therefore N) is not divisible by the divisors. 
We can also use Euclid’s Algorithm or its modifications to test 
divisibility of a number by several primes at once [7,9]. Euclid’s 
Algorithm provides the greatest common divisor of two given 
numbers. If one of these is the number N to be tested and the 
other is a product of some prime numbers, the greatest common 
divisor will be divisible by any of these primes that divide N. 
This method may be used with a large product of primes, say 
up to N!”. The intermediate remainders of Euclid’s Algorithm 
collapse rapidly, but the initial division is an extreme one. We 


Table 1 
Some Divisibility Tests (N = 114043) [3] 


Divisors Reduction Reason 


ig — Of8 = 71 7x il x 13 = 
114 + 043 = 157 31 x 27 =e 
114 — 2(043) = 28 23 x 29 x 3 = 2001 
114 + 4(043) = 286 31 x 43 x 3 = eee 


Der +9448) = 1312 , 19 x 21 = 3 
1140 + 8(43) = 1484 1] x 44 = 7 
1140 + 16(43) = 1828 41 x 39 = 1599 
1140 — 2(43) = 1054 67 x 3 = 201 
114 — 843) = 796 89 x 9 = 801 
bigQe= 9143) = 753 53 x 17 = 901 
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will consider here a less ambitious example to demonstrate its 
possible use in mental calculation. 

First, examine those primes that have been treated so far in this 
section. Working up through the two-digit primes, we have 
derived perhaps the most efficient divisibility tests for the numbers 
aS, 7y 11, 13; 17, 19).23, 29, 31, 37, 41, 43, amdechose of she 
form (10m + 1) and (10m — 1) above this. Therefore, let us 
consider the product 141987 of the primes 47, 53, and 57. This 
number is located near enough to 142000 to ease our calculations. 
We will test the number N = 760603 for divisibility by these 
primes using what we called the modified Euclidean Algorithm. 
We discard easy small prime factors we know are not divisors, such 
as 2 and 5, from convenient positive or negative remainders in 
each step. 


760603 = 5 © 142000 + 50603 
= 5 © 141987 + 4 © 12667 
141987 = 11 © 12667 + 2 @ (5)? ¢ 53 
12667 = 239 e 53 +0 
Therefore, 760603 is divisible by 53, but not by 47 or 57. The 


advantage over long division by each prime increases with the size 
of N and the number of primes included in the product, within 
the capabilities we possess. 

Of course, N can be smaller than the product of the primes. For 


N = 26269 we have 
141987 = 5 © 26269 + 2 e 5321 
26269-= 5 © 5321 —(2)* © 21 


We can stop here since 21 is less than 47, and we conclude that 
26269 is not divisible by 47, 53, or 57. 

The general procedures of nines tests and elevens tests can be 
adapted to other bases. It is occasionally useful, though | have not 
seen it presented or practiced by anyone else, to check a hexadeci- 
mal equivalent by casting out 17’s. For the hexadecimal base, 
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base 16, this is exactly equivalent to casting out elevens in the 
decimal base: 


F29D > (D + 2) -(9 + F) 5 15-24 5-9 3 8 
62109 + 6121109 + 9109 > -9 > 8 


This leads naturally to the question of possible benefits of 
changing bases [9]. Whenever two divisors can be expressed as 
(b + 1) and (b— 1), we can convert the given number to the base 
b and cast out nines and elevens. For example, to determine if 


N = 13949 is divisible by 29 or 31, we can convert to base 30: 
13949/30 = 464 remainder 29 
464/30 = 15 remainder 14 
15/30 =O remainder 15 


Then (15 + 14 + 29) = 58, revealing a factor of 29, and (15 + 
29 — 14) = 30, eliminating the factor 31. 

Converting bases can be awkward when they are not multiples 
of powers of 10. However, we can choose a convenient base that, 
when increased or decreased by one, is a multiple of the divisors 
we are testing. Taylor [9] gives base 50 as such a base for divisors 
3, 7, and 17, as 49 = (7)? and 51 = 17 @ 3. 

Other bases can provide for four divisors less than 100 at once, 
including 300 (since 13 and 23 divide 299 and 7 and 43 divide 
301), 900 (for 29, 31, 17, and 53) and as we have seen before, 
1000 (for 37, 7, 11, and 13). Taylor gives a table like Table 2 of 
convenient bases for divisors less than 100, including all primes 
except 73, 83, and 97. The prime numbers 3 and 9 actually may 
appear in a number of cases here, but are shown in their most 
useful base of 10. 

Actually, we did something reminiscent of this in our earlier 
method for d ending in a 1, 3, 7, or 9. If we regroup the net 
operations in, say, our earlier example of N = 304 and d = 19, 


we find 
301+002+4¢(2)?=19 
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Table 2 
Convenient Bases for Divisibility Tests 


Nines Test Factors Elevens Test Factors 


and for N = 51243 and d = 31, 
S5el-le3+2¢e9-40¢27+3 8 81 = 155 
Lle1-50e3+5¢9= 3] 

Factorization 


Finding factors of integers is often useful in simplifying, for 
example, multiplications of medium-sized numbers. However, for 
large numbers it is surprising in its difficulty; the relative problem 
of factoring large integers compared to multiplying them is 
exploited in certain forms of cryptography nicely described by 
Riesel [7] and Eynden [22]. 

The fundamental theorem of arithmetic states that any factor- 
able, or composite, integer can be represented in one and only one 
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way (barring ones and simple re-orderings) as a product of positive 
prime numbers. Recall that a prime number, or prime, is a positive 
integer other than | that is divisible only by itself and 1. These 
include the numbers 2, 3, 5, 7, 11, 13, 17, and so forth. In this 
section and the next we will examine processes that find factors 
of integers and that lend themselves to mental solution. In fact, 
we will explore this area well beyond its useful scope because this 
represents a fascinating topic for mental recreation and challenge. 

To begin, we can apply our earlier divisibility tests for small 
primes. On the other hand, “Even though a random number n 
usually has small factors (since n is divisible by 2 with probability 
1/2, by 3 with probability 1/3, by 5 with probability 1/5, etc.) it 
is very unusual for n to have only small factors.” [6] 

Riesel [7] remarks that large integers very rarely have many 
prime factors, despite popular perceptions. Moreover, the largest 
prime factor of N would have approximately 62% of the number 
of digits of N. As a cautionary note, the integers in the range we 
work in (or particular ones chosen, of course) may have signifi- 
cantly different characteristics from these. 

To continue, what techniques can we put to bear other than 
divisibility tests?’ One convenient and interesting method (the 
first systematic one, in fact) was developed and used by the 
mathematician Pierre de Fermat (1601-1665) and will find 
the largest factor (prime or not) of N, obviously not greater than 
N!/?. It also has the capability, amazingly enough, of eliminating 
division entirely [6,7,9,10,23—25]. 

First, assuming our number N is odd (the non-trivial case), we 
can locate the midpoint of any two factors, which will both be 
odd, and call this integer x. Then for y < x, 


(x — y) ® Gey) © IN 
ee | 
or 
N= (3) 


Therefore, we must find integers x and y satisfying Equation 3. 
We first set x equal to the lowest value it can be, i.e., the nearest 
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integer above N!”, and increment x in a series of steps. Since 
[((x + 1)? — N] = (x? —N) + (2x + 1), we increase each (x? — N) 
value by (2x + 1), or by the previous increase plus 2. When we 
find that (x* — N) is the perfect square of an integer y, we extract 
(x — y) and (x + y) as factors of N. Obviously, the nearer the 
average of the factors is to N!”, or correspondingly the nearer the 
factors are to each other, the closer the successful x will be to xp. 
In fact, to find the factor a = kN!”, where 0 < k < 1, we require 
(1 — k)N!/2/2k steps, which can be large for small k. 
To illustrate this algorithm, we can take N = 1403: 


a NN 

38 41 

39 4t* 77 = Lie 
40 197 

41 278 

42 361 


Since 361 = (19), then 1403 = 23 © 61, because 42 — 19 = 23 
and 42 + 19 = 61. 

The task of determining whether a perfect square has emerged 
is made much easier by realizing that such a number must have 
its last two digits as 00, el, e4, 25, d6, or e9, where e is an even 
digit and d is an odd digit. In the example above, only 41 and 361 
meet this criterion. 

Actually, we can extend these rules. A square ending in 25 can 
only end in 125, 225, or 625. Also, numbers ending in el and e9 
can only have an even thousands digit if 4 divides e. Endings of 
this sort are: 


e01 e09 
d21 d29 
e41 e49 
d61 d69 


e81 e89 
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Implicit in this, since divisibility by 4 is determined by the last two 
digits only, is the fact that a number of the form (4k — 1) can 
never be a square; squares can only be of the form 4k or 
(8k + 1) [10]. 

These enhancements do not help us in distinguishing 41 and 
361 as squares if we didn’t know. 

In general, we can easily construct a table of values of x? mod m 
(called the quadratic residues of m) which can help; these are 
given for m = 3 to 9 as the third column in Table 3. For example, 
x” mod 9 must leave a residue of 0, 1, 4, or 7, a nines test that 
is passed here by 361 but not 41. Riesel notes that x’ mod 16 gives 
O, 1, 4, or 9 as residues, a test almost twice as efficient. While it 
is very nice for computers to work in this base, it is also a nice 
rule for us. Both numbers pass this sieve as well! 

There is another means of simplifying the entire process using 
modular arithmetic. Since 1403 = 2 mod 3 and for any x, x* mod 
3 = O or 1, then (x* - N) mod 3 = 1 or 2, respectively. Again, 
for any y, y> mod 3 = 0 or 1, so for (x? —- N) = y’, then x’ mod 
3 = 0, giving x mod 3 = 0. For N = 2 mod 3, then, we need only 
consider values of x that are multiples of 3 by the use of this sieve. 
From x = 39 we can therefore jump directly to x = 42, increasing 
(x? — N) by 3(39 + 42) = 243. 

We can construct a table of possible values of x mod m for 
various m. The completed table for some representative values of 
m is given in Table 3. 

In this table, even values of N mod m do not occur for even 
m because multiples of an even number are even and N is odd. 
However, multiples of odd numbers can be even, so both even and 
odd values of N mod m can occur for odd m. 

Looking at Table 3, we note that entries for N mod 6 = 3 and 
N mod 9 = 3,6 do not exist since N mod 3 = 0 in this case. Also, 
we observe that the cases m = 3 and 6 are redundant to the case 
m = 9, and the case m = 4 is redundant to that of m = 8. In 
addition, the cases m = 8 and 9 have relatively few x? mod m 
values for their size, and the results for m = 9 are easily remem- 
bered. In fact, the sums of the possible x mod m values in every 
row of every m, as well as the sums of (x? — N) mod m values, 
are intriguing in themselves. 

In summary, perhaps we can best approach the solution to a 
factoring problem as follows: 
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Table 3 
Modular Sieves for Fermat’s Method of Factoring 


- N) Possible 
N mod m - m x mod m 
0,1,2 0,1 


2,0,0 1,2 (not a 
mult. of 3) 


0 (mult. of 3) 


Ty242 


$,0,3,0 
24,2 


1,3 (odd) 
0,2, (even) 


0,1,4,4,1 4,0,3,3,0 | 0,1,4 


3,4,2,2,4 
#ox1,1,3 
, 240;0,2 


1,4 
oe 
0,2,3 


0,1,4,3,4,1 5,0,3,2,3,0 | 1,2,4,5 (not a 


mult. of 3) 
0,3 (mult. of 3) 


#,2,9;4,9,2 


6,0,3,1, 
1,3,0 


5,6,2,0, 
0,2,6 


7,931,6, 
691.5 


3,4,0,5, 
5,0,4 


2,3,6,4, 
4,6,3 


i220), 


9 


1,3,4,6 


2,3,4,3 


0,2,5 


1 25556 


0,3,4 


0,1,6 


’ 


(table continued on next page) 
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Table 3 
Continued 
— N) Possible 
a m x mod m 


1,3,5,7 (odd) 


2,6 (two + 
mult. of 4) 


1,3,5,7 (odd) 


0,4 (mult. of 4) 


1,8 


0,3,6 (mult. 
of 3) 


2.7 


5,551 2,8 
6005, 

3, 3553036 
4,5,8,4, 
2,2,4,8,9 
202, 

0,0,2,653 
B27 Dyly 
G,8, 1,592 


0,3,6 (mult. 
of 3) 


4,2 


0,3,6 (mult. 
of 3) 


1. Check for divisibility by 2, 3, 5, 7, and 11 using standard 
divisibility tests. 

2. Find xo as the nearest integer larger than N!””. Estimate the 
range required of x; since 11 < (x — y) < 13, a conveniently 
calculated upper limit on x is (N/25) + 6. 

3. Find N mod 9 (by the nines test) to reduce possibilities of 
x. Find N mod 8 to further reduce them. Optionally, find 
N mod 5 and N mod 7. 
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4. Try possible x’s greater than xo, checking the two-digit 
endings of (x? — N) for possible squares. For cases that pass 
this sieve as well, determine if y = (x? — N)! is an integer. 
If so, (x — y) © (x + y) =N. 


This is actually much easier and more powerful than it appears. 
In our earlier example of N = 1403, which appears trivial now, N 
mod 9 = 8 implies that x mod 3 = O. Also, N mod 8 = 3 
(remember to look only at the last three digits) and therefore x 
mod 4 = 2. The latter two conditions, incidentally, amount to 
requiring x to be an odd multiple of 6. In the total range of x from 
38 to 62, only x = 42 and x = 54 are possible, and [(52)* — N] 
does not clear the two-digit ending sieve. The result x = 42, 
y = 19 is therefore easily obtained. 

If we had memorized or reconstructed the rules for m = 5, we 
could have found that x mod 5 = 2 or 3. Then x = 42 is the only 
candidate, the next possible value, 78, being beyond the range of 
x. For m = 7, if we bothered to recall it, we would find x mod 
7 = 0, 2, or 5, and the next possible value after x = 42 would 
be 114. 

As mentioned, Fermat’s “difference of squares” method of 
factoring works from y small toward y large, so two factors a and 
b far from each other will take a long time to locate. If it seems 
that we are not having any luck, we can try multiplying N by a 
factor k. For odd k, we will find the factors a and kb quickly if 
a is approximately equal to kb. The case of even k’s is a little 
complicated, for then kN will be even. For k = 2, KN cannot have 
a midpoint of two factors. If we multiply x and y by 2, then 
(2x — 2y) © (2x + 2y) is a valid factorization, giving k = 4, but 
the factors are the same relative distance apart. For k = 8, we then 
have the valid factorization (4x -— 4y) @ (2x + 2y), giving a fast 
solution for a approximately equal to 2b. In general, any even 
value of k should be a multiple of 8 when a is an even multiple 
of b. A case where a/b is approximately 3/2, of course, would 
benefit from setting k = 6, the case a/b approximately 5/3 by 
letting k = 15, and so forth. 

Given that we generally have no feel for the relative sizes of the 
factors, we can abandon the factoring process after a while and try 
multiplying by a new factor. Riesel [7] suggests multiplying by a 
composite k containing various small divisors (such as 24). There 
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is a legendary feat of factorization by Fermat, often quoted as 
requiring a primality test lost in history, wherein the number 
N = 100895598169 was reduced to 112303 x 898423. It turns out 
that 8N = (898424)* — 898424 = 898424(898424 - 1) = 8 x 
112303 x 898423. 

There are situations where we can eliminate most values of x 
in Fermat’s method without constructing or recalling the modular 
sieves we have developed. While this falls very near the domain 
of properties of particular numbers, it is of such great potential 
that I feel it is worth discussing. 

It is known [7] from A. M. Legendre (1752-1833) that odd 
integers N of the form N = a" + b" for any integer n and of the 
form N = a® — b" for any odd n, where gcd(a,b) = 1, have prime 
factors of the form p = 2kn + 1, excluding prime divisors of 
algebraic factors of N of the form (a™ + b™) for m < n. Of course, 
for N = a" — b" for n even, we already have factors as (a" + b"*). 

We will discuss the clause concerning algebraic factors in a 
moment. In general, this sieve may appear to limit possible values 
of x in Fermat’s method to 1/n of their number otherwise, but it 
actually limits them to 1/(2n7), a significant advantage. In fact, we 
need only consider values of x for which 


+ | 
mod (2n7) 


The occasional number of the given forms (a + b") benefit greatly 
from this sieve. There are more of them around than you might 
expect, particularly since b can be as small as 1. 

As an example, consider N = 1027 = (10)° + (3)°. Then x = 
514 mod 18 = 10 mod 18. Since our values of x would normally 
begin at 33, we can now jump directly to x = 46 and begin 
incrementing by 18: 


x xa 
46 1089 


We find 1089 = (33)? immediately, so 1027 = 13 @ 79. 
Now we return to the matter of algebraic factors. These 
are factors of algebraic expressions, such as 3xy + 2y — 3x —2 = 
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(3x + 2) © (y—1). Of course, if N can be expressed as a factorable 
polynomial, we immediately have two factors. Now, since 


a2 — b? = (a + b) @ (a—b) 

a’ — b? = (a—b) @ (a” + ab + b’) 

ee bte (e + bth Ore? ab?) 

=a 7p) ©" (a + b) ® (a — b) 

a’ — b? = (a—b) @ (at + a’b + a’b? + ab? + bt) 
and 

a’ + b? is unfactorable 

a? + b? = (a + b) © (a2 — ab + b’) 

a* + b4 is unfactorable 

a? + b> = (a + b) @ (at — ab — a*b? — ab? + bt) 


then most of the expressions (a + b") or (a" — b") reveal two 
factors immediately. The factor (a + b) is associated with the 
expression (a" + b"), and the factor (a — b) with the expression 
(a" — b"). Therefore, we could instantly have known earlier that 
one of the factors of 1027 = (10)? + (3)? is (10 + 3). 

Legendre’s theorem is powerful because it defines a sieve for all 
prime factors other than those dividing these algebraic factors of 
the form (a" + b"), and as we discussed, most numbers contain a 
variety of prime factors. Riesel cautions us to test if prime factors 
derived from algebraic factors are multiple factors before we search 
for those of the form (2kn + 1). 

Let us examine a trivial case, N = 1001 = (10)’ + 1. We 
immediately know that (10 + 1) = 11 is one factor, and it 
accordingly is not of the form (6k + 1). The only algebraic factor 
of this type is (a + b), and no other prime divides it. Dividing N 
by 11, we find the cofactor N' = 91 (which is not divisible by 
11 again) and we know its prime factors now to be of the form 
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(6k + 1). We can test k = 1, 2, 3, etc., finding 7 and 13 or, more 
generally, we can use Fermat’s method with x = 46 mod 18 = 10 
mod 18. Then we find that (10)? — 91 = (3), giving us the 
remaining prime factors 7 and 13. 

Considering a less trivial case, let N = 2581 = (50)? + (9). 
There are no algebraic factors here, so all prime factors are of the 
form (4k + 1). Using Fermat’s method, we know then that 
x = 1291 mod 8 =3 mod 8. We would normally begin this method 
with x = 51, and since 51 mod 8 = 3, we persist, increasing x in 
each step by 8: 


x x? —N 
51 20 
59 900 = (30) 


Therefore, 2581 = (59 — 30) @ (59 + 30) = 29 © 69 Tf as 
recommended earlier we had memorized the modular sieve from 
Table 3 for m = 9, then the result N mod 9 = 7 would imply that 
x mod 9 = 4 or 5, catapulting us immediately to x = 59 as our 
initial value. 

There is a variation of Fermat’s method called, by Vaes, the method 
_of remainders [10]. Setting y = (N - 1)/2 and x = (N + 1)/2, 
giving N = x’ — y’, then a remainder r obtained by dividing y by 
p will give a remainder (r + 1) when x is divided by p. Therefore, 
if we divide y by p and find r, and (2r + 1) is divisible by p, then 
N is divisible by p. 

For example, for our earlier N = 1403, we divide (N — 1) by 
2 (an immediate advantage) to arrive at y = 701. Then we begin 
with Equation | to rewrite the value of y in each step, producing 
various remainders r. We have 


y= (26) 25 
= 25 o2/ + 2 
= 24 ¢ 28 + 29 
= 23 © 29 + 34 


a 


Primitives 57 


and we continue, where the remainder r is the last number in each 
step, until (2r + 1) is a multiple of one of the multiplied terms 
in that step. In the last step shown, 2 © 34 + 1 = 69, a multiple 
of 23, so 23 is a factor of 1403. Perhaps finding whether a number 
is divisible by another is easier than finding if it is a square, but 
the memory requirements are greater at each step. We can ease 
the task if, instead of adding the squares of 1, 2, 3, and so forth 
to our original remainder in each step, we simply add the odd 
numbers 1, 3, 5, 7 to the remainder in the previous step or, better 
yet, add twice the odd integers cumulatively to our original 
(2r9 + 1) value: 


to = 25 

2t— + 1 = 51 

2r, + 1 = 51 +201 = 53 
277+ 1=53+2°3 =59 
27, + 1 =59+2°5=69 
etic. 


In a different vein, we can use triangular numbers instead of 
squares in Fermat’s algorithm [9,10]. Triangular numbers represent 
the total number of items in the rows of a triangular array up 
through the tth row, such as in the pin arrangement in bowling. 
The triangular numbers are given, then, as 1, 3, 6, 10, 15, 21, ete. 
and are given by the expression 


ux + | 
(oO a 7, 8,... 


Now, since 


ae Py yy > 1) Ge = ye (e+ y tt) 
2 2 2 


OF 
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then if we find a triangular number t, that after subtraction of N 
results in another triangular number t,, we have factored N. The 
denominator 2 in the factored expression for N divides whichever 
term in the numerator comes out even. 

As we did with squares, we begin with the nearest t, above N; 
x is found here by extracting the square root of 2N. Since 


x(x + 1) (ke - D&-1+D) 
, a 2 


xX 


Or 
ty = bey 


we can increase the previous value (t,_; — N) by x in each step. 
Of course we knew that, since we form consecutive triangular 
numbers by adding another row in the bowling pin arrangement 
and the xth row contains x pins. Therefore, we do not have to 
keep track of our value of x as we proceed. In our earlier example 


of N = 1403, we have 
(2806)! is approximately 53 


(53 © S42 = ee 


and we collapse immediately as 28 = (7 © 8)/2 = t7. The factors 
of 1403 are then (53 — 7)/2 = 23 and (53 + 7 + 1) = 61. We seem 
to bear out Taylor’s claim that the great advantage of using 
triangular numbers over squares lies in the reduction in the 
number of steps involved. Experimenting with various numbers 
reveals the startling truth of this for factors not exceedingly close 
to one another. 

Let’s try factoring another number to show the values added in 
each step: 
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N = 21449 
(2N)!” is approximately 207 
207° 208 200 © 215 + 56 


= 21528 
2 4 
x (or adder) t, -N 
207 79 
208 79 + 208 = 287 
209 287 + 209 = 496 


but 496 = (31 © 32)/2. Then the factors are (209 — 31)/2 = 89 
and (209 + 31 + 1) = 241, showing an amazing convergence to 
a solution with such widely differing factors. 

How do we determine if a number is triangular? Triangular 
numbers end in the digits 0, 1, 03, 53, 5, 6, 28, or 78, and the 
first such number in the last example is 496. The single-digit 
endings listed here can be preceded by any digit, so we unfortu- 
nately end up with 44 possible two-digit endings instead of the 22 
we had for squares. Now t, mod 9 = 0, 1, 3, or 6, a test easily 
performed by casting out nines. Since 496 passes this sieve as well, 
we can explicitly check it by finding (2 ¢ 496)!” as approximately 
31. Then we find that 2(496) = (31)? + 31, or equivalently that 
2 © 496 = 31 © 32, and we are done. 

We can construct a table, as for Fermat’s method, of possible 
values of x mod m for various moduli m. This is given in Table 4 
for various moduli. 

Again, even values of N mod m will not occur for even m, as 
N is odd and multiples of an even number are even. We also 
eliminated N mod 6 = 3 and N mod 9 = 3 and 6, as these would 
then be divisible by 3, a prime factor considered to have been 
tried already. 

Overall, we are not as fortunate here as we were for Fermat's 
method. The cases m = 3 and 6 are again redundant to the case 
m = 9, and perhaps in general the most efficient procedure is to 
simply memorize the case m = 9, which is not difficult. In our 
previous example, 21449 = 2 mod 9, so x = 2 or 6 mod 9. The 
first such x > 207 is 209, and the next occurs at 213. To jump from 
209 to 213 we would have added 4 © 211.5 = 846 to arrive at 1342. 
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Table 4 
Modular Sieves for the Triangular 
Number Method of Factoring 


(t, — N) Possible 
t, mod m i m x mod m 


4 


0,1,3 


1,2,3 (not a 
mult. of 4) 


0,1, 5,108 4,0,2,0,4 


3,4,1,4,3 
Ly IiOnoee 
Ls 2iAgel 


031,350,490 5,0,235,5,5 


1,2 445 


0,2,3,3 


6,072, 5, 
2,0,6 
5,6,1,4, 
Loe 
459058, 
0,5,4 
Sehta Bad 
6,4,3 
Zeal, 
55 Sind 
1,2,4,0, 
#2, 1 


0,1,5,6 


1,2;4,5 


2,354 


0,2,4,6 
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Table 4 
Continued 
Possible 
N mod m x mod m 


d Dees 
1,054,¢ 


5,6,0,3, 


’ 


’ 


3,4,6,1, 
5,2,0,3 


’ 


4,5,6,7,8 


6,2,6,0,5 
4,5,/,1, 
5,1,7,5,4 
pMaRonen 
5,8,9,5,2 
1,2,4,7, 
2,/,4,2,1 


There is another method of factorization for numbers of the 
form (4k + 1), where k is a positive integer. It has been proven 
that a prime number of this form (i.e., 1 mod 4) can be expressed 
as the sum of squares (a” + b’), with gcd(a,b) = 1, in one and only 
one way; otherwise, the number is composite. Fermat first dis- 
covered this relationship which was proved a hundred years later 
by Euler. Fermat called this fact “the fundamental theorem on 
right-angled triangles,” as for a prime p there is only one such 
figure with integral sides and a hypotenuse equal to p'/*. The 
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method of factoring when there are two sums of squares is termed 
Euler’s method, although it was given much earlier by Frénicle and 
Mersenne [23]. 

] mention this technique largely because of its use by at least 
two lightning calculators in this century [3]. In my opinion, its 
main disadvantages are that it can be used on only half of all odd 
numbers, it requires checking every possible sum of squares (albeit 
sieved), and it requires subsequent extraction of the factors from 
the squares if the number is composite and has at least two such sums. 

There are two basic methods for finding pairs of squares that add 
to a number N, which we term for two sets of squares 


N=a@@+b=c2 +d? 


The first method parallels Fermat’s method of factoring, where 
we now find whether (N — x’) is a square rather than our previous 
quantity (x? — N). We can test every value a between N!” and 
(N/2)"2, or 0.7N!, at which point we begin duplication with x 
and y interchanging roles. In analogy with Fermat’s method, we 
begin with x just less than N!/” and decrease x as we proceed. We 
need only add (2x + 1) in each step, or simply increase the adder 
in the previous step by 2: 


N —- x? = N-—-(x + 1)? + 2x + 1 


The two-digit endings for squares can be used with advantage here. 
Let us take the example N = 1433, which is of the type allowed 
since 33 = 1 mod 4. 


zm. We x? 

B/| 64 

36 137 = 64 + 73 
5 203 =15/ + 71 
34 Lt 

33 344 

By: 409 

6 I 472 

30 533 

29 592 

28 649 

ae 104 


26 OW 
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We stop because 757 is larger than (N/2). We see as we mentally 
trace these steps that there is one and only one sum of two 
relatively prime (gcd(a,b) = 1) squares that add to 1433, (37) + 
(8)? = 1433, revealing that 1433 is prime. 

Let’s try N = 1457: 


a. Niwa 
38 13 
37 88 
36 161 
35 232 
34 301 
33 368 
52 433 
31 496 
30 557 
29 616 
28 673 
27 728 


and we find that 1457 is composite, but we have no sum of squares 
to extract factors. Actually, this situation will occur whenever any 
prime factor of N is of the form (4n + 3) and is of an odd power; 
here, 1457 = 31 © 47 and both factors are of this form [26]. 
Fortunately, this is not always the case. For N = 1537, we have 


Lz N= 
39 16 
38 93 
37 168 
36 241 
35 312 
34 381 
38 448 
32 513 
31 576 
30 637 
29 696 
28 753 


27 808 
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We find that 1537 = (39)? + (4)? = (31)* + (24). Now we wish 
to determine factors of 1537 from this information. 

The basic procedure for N = a? + b? = c? + d? = (possibly other 
sums of pairs of squares) is to first find the following values f and 


g [10,23,24]: 
f=—a4—-c 
g=d-b 


given here by f = 8 and g = 20. Now we find the quantity (2ef + 
f?) = (2bg + g’) = 560, and we divide this by fg = 160, giving 
h = 3.5 in this case. Then 


N = (f° + g’) @ (1 + h’)/4 
= (464 © 13.25)/4 
= 29 e 53 
We could have directly used the relation, 


_ {fa - cP + (b - d)’] © [a +c)? + (b- 


N 
4(b — d)’ 
divvying up the denominator into products that divide the terms 
in the numerator. 

Of course, we can use modular arithmetic sieves to decrease the 
number of possible values of x. Construction of a table such as 
Table 3 for the quantity (N — x”) would show that for N = 1537 = 
7 mod 9, then x = 0, 3, 4, 5, or 6 mod 9, small relief in this case. 
However, for an N = (5k' + 2), as here, x and y equal 1 or 4 mod 
5, or in other words are of the form (5p + 1). The fact that N 
is of the form (4k + 1) does us no good, since all x mod 4 are 
possible, and since N = (4k + 3) eliminates all x mod 4 we begin 
to see why we are limited to the former type. 

We can also utilize two-digit endings of squares to determine 
possible sums of two squares that will give the last two digits of 
N. An analogous situation will be presented later, where this 
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sieve is used in finding possible differences of two squares in 
Fermat’s method. 

As a comparison for this example, Fermat’s method without the 
use of sieves takes only two steps before the factors 29 and 53 are 
found, although for widely differing factors or for prime numbers 
the methods are not as dissimilar. 

For those of us still bothered by our inability to factor N = 1457 
earlier with this method, despite showing it to be composite, we 
can turn to a more general form of Euler’s factorization [7]. Here 
we find two representations of a number N in the form (a* + Db’) 
and (c* + Dd’). Since all numbers cannot be expressed in this 
form, we might veer off into this method when we find an entry 
(N — x’) which is a multiple of a square greater than 1. In our 
earlier sequence for N = 1457, we find in the second step that 
N — (37) = 88, so that N = (37)? + 22(2)2. Then we start a new 
series of tests to find a square result of (N — 22x’), beginning with 
x as the nearest integer less than or equal to (N/22)!. Since we 
are decreasing in each step the initial value of x = 8, we have at 
most eight steps. In each step we can add 22(2x + 1) to our 
previous value of (N — 22x’). We begin with 


a N — 22x? 
8 49 


and we stop immediately, having that 1457 = (37)? + 22(2)° = 
(7)? + 22(8). 

To find the factors, we use the more general form of our previous 
technique—we find gcd(N,ad — bc) or gcd(N,ad + be), each 
representing a factor. Here (ad — be) = 282 and (ad + be) = 310, 
which offers, after a quick solution using Euclid’s Algorithm 
and its variations, the factors 47 and 31. (Try Fermat’s method on 
this number!) 

As a second general method, we can find pairs of squares that 
add to N of the form (4k + 1) by finding pairs of triangular 
numbers that add to N' = (N — 1)/4 [10]. This derives from the 
algebraic equality, 


tbe ire t-weeg PEA 
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Considering our earlier example of N = 1537, we have N' = 384. 
Then we find (N' — t,) for x up to N'!”. 


N'-t 
384 
383 = 384-1 
381 = 363 -— 2 
378 = 381 -—3 
374 
369 
363 
356 
348 
339 
329 
318 
306 
293 
219 
264 
248 
238i 
Dad 
194 


a 
ODMDINUDNMNAWNHNHOWODOQVANABWNHKOI* 


From endings and/or modular arithmetic sieves, we find that 
378 = to7 and 231 = t»,, again giving us the intermediate result 
1537 = (39)? + (4)2 = (31)? + (24). 

There is an intriguing if not foolproof test for determining 
whether a given number N is prime if N is of the form (4k + 1) 
and ends in 3 or 7 (e.g., 1433) [9,10]. Sadly enough, even when 
the trivial case of N ending in 5 is eliminated, this test applies to 
only one in four odd numbers. 

First, we subtract the nearest square ending in 5 from 2N, 
leaving a remainder rp. These squares actually must end, as we 
know, in 25 and are the squares of numbers ending in 5, or in 
other words (5a)’ with a odd, so our shortcut for squaring 
such numbers helps a great deal here. Then we find the follow- 
ing quantities: 
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ry =) + 100(a— 1) 
r =r, + 100(a - 3) 
r; =r + 100(a — 5) 
etc. 


Now if one and only one square appears among _ these 
values, N is divisible by this square or is prime. Otherwise, N 
is composite. 

For N = 1433, we have 2N = 2866. The nearest square ending 
in 25 is found by determining a value b such that b(b + 1) is 
nearest without exceeding 28. Therefore, b = 4 and (45)? = 2025. 
Continuing, 


a=9 

fo = 2866 — 2025 = 841 = (29)? 
r, = 841 + 800 = 1641 

r. = 1641 + 600 = 2241 

r; = 2241 + 400 = 2641 

ry = 2641 + 200 = 2841 


There is only one square here, and 841 does not divide 1433, 
so 1433 is very probably prime according to this method. In fact, 
1433 is prime, as we have seen. 

For three-digit numbers this is very fast: 


N = 713 

2N = 1426 

5)" = 1225 
tf = 201 


es 801 
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No value of r is a square, so 713 is composite (713 = 23 @ 31). 

Composite integers are identified without fail in this test, but 
some numbers indicated as prime are actually composite. Taylor [9] 
notes that these latter anomalies “.. . are in general those integers 
which are multiples of primes of one of the forms (20k + 13) and 
(20k + 17) and also of one or more of the squares of 3, 7, 9, 11, 
23, 27,... (ie., (20m + 3,7,9,11)*).” They include 153, 333, 477, 
657, 833, and higher numbers. If the factors 3, 7, 11, and 13 are 
eliminated first, however, an anomaly will not surface until N 
equals 8993. 

In closing, I mention for those interested in pursuing this 
subject on a computer or programmable calculator that summaries 
of factorization methods specifically designed for these instruments 
can be found in Riesel [7] and Blair [27]. Other methods suitable 
for hand calculations (or the above instruments) are given by 


Taylor [9] and Dickson [10]. 


A Factoring Game 


Now I would like to take an extreme case of factorization to 
illustrate additional tests. While this will turn out to be unrealistic 
to perform on the spur of the moment, we can consider this 
factoring of a very large number to be a kind of game, and it’s 
actually not a bad one at that if you’re interested. 

We would like to factor a given large number like 125,869. 
Since N is not divisible by 2, 3, 5, 7, 9, or 11, we begin using 
Fermat’s method (my choice) with x = 355. Chapter 3 provides 
convenient methods for simple square root extraction. Here we 
would take (1258)!/" as roughly 35; then [1258 — (35)?]/(2 @ 35) 
as roughly 0.5. Appending 5 to 35 and checking the decimal point 
location, we arrive at 355 as N!” to the nearest integer. 

Unfortunately, we now find that the upper limit on x is 
(125869/25 + 6) = 5041, leaving an absolutely huge range 
to consider! 

Undaunted, we find N mod 9 by casting out nines and arrive 
at 4, implying from Table 3 that we need consider only values of 
x = 2 or 7 mod 9. We therefore decrease the number of possible 
values of x to roughly 521. We can then find N mod 5 = 4, giving 
x = 0, 2, or 3 mod 5, and reducing the choices of x to around 313. 
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Finally, N mod 7 = 2, giving x = 2, 3, 4, or 5 mod 7 and a net 
of about 179 choices. 

Actually, since the formula for the upper limit on x is roughly 
inversely proportional to the lowest prime p not tested by standard 
divisibility tests, 


we can see that we really suffer in considering all primes in 
Fermat’s method down to p = 13. Let us set p equal to 101 and 
perform two-digit divisibility tests for all primes below this if 
required. Then x,,,, = 679 and our original set reduces from 4686 
to 324 values of x. Then applying our tests N mod 9, 8, 5, and 
7 reduces the number of possible values of x to about 12. 

We find no primes below 100 to divide N. However, to 
sequence through values of x to find the approximately 12 possible 
values from 355 to 679 that pass the sieves, and for each one 
square x, subtract N and check for endings of those that may be 
squares, explicitly finding square roots of those, is tedious still. Can 
we use more detective work to help us out? 

Yes, we can. Let's see if we can find other sieves in a different 
manner. We know that N ends in 69 and that squares end in 00, 
el, e4, 25, d6, or e9, where e is even and d is odd. What possible 
difference of squares end in 69? A quick check of Table 5 produces 
the number of possibilities for differences of squares ending in a 
particular digit. 

Obviously, we could have been luckier here. The four possibil- 
ities for a difference ending in 9 can be easily deduced to be 
(00 — el), (e4 — e5), (25 — d6) and (e9 — 00). 

For numbers ending in 69, we can reduce the choices to (25 — 
56) and (69 — 00). We can now find what x must end in to arrive 
at these endings. It turns out that for each of the 22 possibilities 
of two-digit endings for a square x’, there are four possible two- 
digit endings of x, except for 00 and 25, which have ten (numbers 
ending in 0 and 5, respectively). For the former, when one two- 
digit ending, say x), is found, the other three can be found as 
1100 — x, 1, 150 — x,| and 150 + x,!, where here the vertical 
bars indicate the absolute value of the term inside. This follows 
from the equations, 
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Table 5 
Differences of Squares Ending in Particular Digits 


Number 
of 


Possibilities 


Og 
a. 
— 
aps 


9 
8 
7 
6 
5 
4 
3 
2 
l 
0 


nN F&F NY NY FS DW HS LY LH 


(100 — x,)? = 10000 — 200x, + x? 
(50 + x,)? = 2500 + 100x, + x} 


For example, we know that (13)* = 169, which ends in 69. 
Therefore, other two-digit endings that when squared produce an 
ending of 69 are (100 — 13), (50 + 13), and (50 — 13). 

We can, then, find possible endings of x for each ending of x’. 
We are in reality performing a manual mod 100 sieve. 


x’-ending x-ending 
69 13,39 (05507 
25 05,15,25,35,45;55;65, 75,85}95 


Notice that we can already see that the sieves x mod 8 = 1, 3, 
5, or 7 and x mod 5 = O, 2, or 3 are useless here, since they can 
be determined by these endings. More work is still necessary. 
While it may seem lucky to have isolated the endings to particular 
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numbers instead of e4, for example, the value 25 with its ten 
entries above hurts us considerably. 

In reality we don’t need to construct the above x-endings as a 
table in our mind, but rather find each x-ending one at a time and 
operate the rest of our procedure on them individually. We then 
end up sequencing one line at a time through Table 6. Notice that 
355 <x < 679. 

Now by using endings we not only have reduced the number 
of possible values of x to four, but we have also identified them 
without testing intermediate values. At this point, if we want to 
avoid squaring the remaining x-values, it is appropriate to check 
y’ for any reduction of the y-values as well. 

Fortunately, this is very easy. Since adding N mod m is the same 
thing as subtracting (m — N) mod m, then we can use Table 3 for 
x mod m and use the values derived for the complement of N mod 
m, thereby giving us general results for a number (y? + N) mod 
m. Therefore, without actually doing any divisibility tests we 
immediately know that 


y mod 9 = 0,3,6 from entry for N mod m = 9 — 4 = 5 
y mod 8 = 2,6 from entry for N mod m = 8 — 5 = 3 
y mod 5 = 0,1,4 from entry for N mod m = 5 - 4 = 1 
y mod 7 = 0,3,4 form entry for N mod m = 7 - 2 = 5 


Of these, we are interested in y mod 5 and y mod 8, since 
divisibility by 5 and 4 is obtainable from the two-digit ending. 

Looking at endings for y?, we can deduce that a y’-ending of 
56 arises from y-endings of 16, 34, 66, and 84. Of these, the above 
two sieves limit y-endings to 34 and 66. 

We can now take each of the four values of x remaining after 
our sieve and find possible values of y. Remember that (x — y) > 
100. In addition, we know that (x? — N) = y’, and this must hold 
up to all divisibility tests. 

For x = 605, we know that (x? — N) mod 9 = [(2)’ — 4] mod 
9 = 0; because of the earlier x mod 9 sieve, all of these values of 
x will give this result. Therefore, only y = 234 is possible, since 
the nines test on the ending 34 implies an initial digit of 2 tor the 
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Table 6 
Modular Sieve Results for N = 125869 


Sieve Additional Sieve 
x mod 9 = 2,7 | x mod 7 = 2,3,4,5 


range 0 < y < 505. An ending of 66 produces no initial digit in 
this range. As y = 234 seems a likely possibility, we cast out 
elevens to check it again, producing (x’ — N) mod 11 = [(0)? — 7] 
mod 11 = 5. The y-value 234 is now eliminated because (234)? 
mod 11 = 9; thus x cannot be 605. 

For x = 515, again only y = 234 is possible and is eliminated 
by the elevens test. 

For x = 425, we find y = 234 does indeed pass the elevens test. 

For x = 565, we find that y = 234 passes the elevens test as well. 

We can perform other divisibility tests if we still want to avoid 
squaring three-digit numbers. As mentioned earlier, a good one is 


for 37. 
(x? — N) mod 37 = [(425)? — 125869] mod 37 
= [(-19)? — (014,-019)] mod 37 
= [367 — (-19 + 14)] mod 37 
= 366 mod 37 = 33 
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As a reminder, we dropped multiples of 111 out of each triplet 
and added each triplet to the triplet to the right of it. Now (234) 
mod 37 = (12)? mod 37 = 33 as well. 

If we repeat this for x = 565, we find (x? — N) mod 37 = 31 
and the failure of y to pass this sieve eliminates this possibility. 

We've danced around the result for quite a while now and we 
know that x = 425, y = 234 is strongly indicated, since out of two 
cases like this the odds of one passing the divisibility test for 37 
without being a solution is one in 18.5. We could square this out 
if we wanted, but let’s continue to be stubborn about this and 
assert with great confidence that (425 — 234) = 191 and (425 + 
234) = 659 are factors of 125869, which indeed they are. 

Without using modular arithmetic sieves, how many steps 
would we have had to perform using Fermat’s method to find these 
factors? We earlier had the expression (1 — k)? N!/?/(2k) for the 
number of steps, given that the lower factor a = kN!”. Here 
a = 191 and we find k = 0.538, giving 70 steps (after the initial 
value of y is obtained). Using triangular numbers, which shine 
when factors are disparate, we require 18 steps without sieves. 
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Roots 


The real test of ability to do square, cube or any 
other root is, in my view, to have a number 
proposed that is not an exact power, and to be 
asked to give the answer to several decimals; but 
this type of question you will hardly find in the 
published records. 


A. C. Aitken (1954) [1] 


A\tthough extracting integer roots of powers is a common 
practice of lightning calculators, deriving roots of numbers 
that are not perfect powers is historically much more rare. 
This chapter is concerned with the much more practical 
problem of determining roots of such diverse numbers. None- 
theless, we will first touch on methods for perfect powers, if 
only to satisfy curiosity on the subject. 


Roots of Perfect Powers 


Roots of perfect powers, particularly odd powers, offer 
distinct advantages to the calculator. Briefly, it is well known 
that for a power of order (4k + 1), with k a positive integer, 
the root will contain the same units digit as the power. For 
a power of order (4k + 3), the root will contain a unique units 
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digit for each unit digit of the power. For a three-digit root 
of any odd perfect power, the leading digit can be inferred 
from memorized or estimated ranges of such powers. With the 
knowledge of the units digit, the middle digit can be derived 
by casting out nines and/or elevens. 

To illustrate, suppose we are given 6657793506607 as the 
fifth power of a number and our task is to find the fifth root. 
Now, the given number is about 670 x 10!° and we find (or 
memorize) ranges for the first digit of the fifth root. We see 
that (300)? = 243 x 10!° and (400)° = 1024 x 10!9, so the 
answer is a three digit number with 3 as its first digit. Since 
the order of the root, 5, is of the form (4k + 1), we 
immediately know that the last digit is the same as that for 
the given number, or 7. To determine the middle digit a, we 
use the nines test and, if needed, the elevens test: 


6657793506607 mod 9 = (6+6+5+7+7+9+3+4+5 
+0+6+6+0+7) mod 9 
= 
We know, then, that (3a7)° mod 9 = 4. 


Now we want a number less than nine that, when raised to 
the fifth power, will leave a remainder of 4 in the nines test. 


Nines test result 
Number n on nm 


COMIN nh WN — © 
CORONA 10 NK OC 


We could have memorized this table or generated it on the spot. 
To do this, we remember that we can always reduce intermediate 
values in our calculation by the nines test. For example, for n = 5, 
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n’ mod 9 = (5)? mod 9 = 7 
n* mod 9 = (7)? mod 9 = 4 
n? mod 9 = (5 © 4) mod 9 = 2 


The only ambiguity in this test is for a result of zero, which 
corresponds to n = 0, 3, and 6. The elevens test can then be used 
to distinguish between these cases. For cube roots, the elevens test 
is preferable up front, because it provides unique results for every 
*<\)1. 

In our case, we have now deduced that 3a7 mod 9 = 7, or 
(3 + a + 7) mod 9 = 7. Simple inspection shows that a = 6, so 
we arrive at 367 as the fifth root of 66577993506607. 

For roots containing more (say, m) digits, the logarithm may be 
calculated of the first (m — 2) digits of the power. The result is 
then divided by the order of the power, as we do when using 
logarithms, and then converted back into the initial digits of the 
root. The last two digits are then derived from the uniqueness of 
the units digit and casting out nines or elevens, as before. 

For the practiced calculator, the task can be simplified by 
memorizing tables of two-digit or three-digit endings of powers. 
For most roots, the endings are not unique, though, and require 
a subsequent process of elimination. Nonetheless, the huge ranges 
of powers that can be considered can make the extraction of 
perfect roots seem amazing. For example, extracting a three-digit 
cube root, which as described above is trivial, can be performed 
on numbers up to (999)? = 997002999. All three-digit roots of odd 
powers are equally trivial, and the ranges of numbers become 
incredible for high powers. Smith [2] provides an excellent review 
of these techniques and the use of endings to ease multidigit 
integer roots. 


Particular Square Root Methods 


Since we rarely find ourselves extracting roots of powers known 
to be perfect, a more desirable algorithm would extract noninteger 
roots as well. A few of the better methods of obtaining a 
noninteger root (and to be precise, an integer root as well) were 
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described by Alexander Craig Aitken, considered by Smith to be 
“one of the greatest mental calculators who ever lived [2].” 

In a passage from Aitken’s address to the Society of Engineers 
[1], part of which can be found in Smith [2] as well, we find a 
presentation of six ways an experienced calculator might approxi- 


mate the square root of 51 (or 7.141428428543 ...): 


1. Take an initial estimate of the square root of 51, or (51), 
as 7. The average of 7 (or 49/7) and 51/7 provides a much 
better estimate, 50/7 = 7.1429.... 

2. Take 50/7 as the new initial estimate. The average of this 
and 51/(50/7) is 7.141429. ... 

3. Take the repeating decimal 7.141414...as a new initial 
estimate. Referring to our discussion in Chapter 2, we know 
this to be the fraction 707/99. The average of this and 
51/(707/99) is 7.14142842857 ..., accurate to eleven digits. 

4. Return to our estimate of 50/7, which was the average of our 
initial estimate of 7 (or 49/7) and 51/7. Their ratio is 49/51 
and each “deviates” by one part in 50. Reduce 50/7 by one 
part in 2(50)?, or one part in 5000, arriving at 4999/700 = 
7.1414285. ... 

5. Divide the interval from 49 to 51 into fourths and 
multiply 7 by the third quarter divided by the first, 
giving 7(504/494) = 7.141414141414.... Alternatively, 
find (51/7) © (494/504) = 7.141442715700. ... Their aver- 
age is even better, being, in fact, that found in method 3. 

6. Using “subtler and more powerful approximations still,” 
reduce 50/7 as in method 4 by one part in 4999% rather than 
one part in 5000, giving 7.141428428557 ..., “so commit- 
ting an error of 1 in 500,000,000,000. This is an extreme 
approximation for [a] square root; and I have never gone 
beyond it in mental calculation.” 


While Aitken does not elaborate in his presentation on the 
bases for his methods, we can with some effort reconstruct the 
reasoning behind them. To begin, we make use of a concept 
formulated by Heron of Alexandria (c. 50?) and most likely by the 
Babylonians of 1700 B.C. [3,4]. Aitken (and Smith) did explicitly 
refer to this method. 


se 


; 
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We reason here that if a number N is deduced to have an 
approximate root a (say, one that is too low), then the quantity 
N/a provides a number that will lie on the other (high) side of 
the root N!?, Therefore, a better approximation to N!/? will be 
the average of a and N/a. This procedure can be iterated and is 
of second order, i.e., the number of correct digits approximately 
doubles with each iteration [5-7]. We then have the iterative 
relation with n = 0,1,2,... 


a, + (N/a,) 
ane} = ————— (4) 
Z 
N’ —-af 
= a, + ay 
2a? 


This is actually a special case of the well-known Newton- 
Raphson method of finding simple, real roots of a real equation 


of the form f(a) = 0 [3,8—11]: 
f(a,) 
~ fi(aq) 


where f'(a,,) denotes the first derivative of f(a) with respect to a 
evaluated at a,. For f(a) = a? - N = 0, 


(5) 


an+1 = An 


an+1= Ap + an (6) 


which reduces to Equation 4 when the order of the root p = 2. 
Aitken derived his first three estimates from his initial approxima- 
tions using one iteration of the Newton-Raphson method. 

There are several ways of deriving the Newton-Raphson rela- 
tion (Equation 5). A nongraphical derivation follows from the 
expansion of a function f(a) into a Taylor series: 


fla) = flay) + F(a,) * (aay) +" » 


(a — a,) 
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We can truncate this after the second derivative term and set 
f(a) = 0. Casting the formula in an iterative form to approach f(a), 
we then have 


f(a; ) 


° (ans = ae =@ 


f(a) + flag) ® Gaeta) + 


which can be rewritten as 


bee -f(a,) 
flay) + (eqs — ah)? f'faq)/2 


This reduces to Equation 5 when the second derivative term in 
the denominator is neglected. 

Certainly more satisfying (if less general) is a simpler derivation 
for roots of order p. Here, for an error e,, a root a = a, + e,, and 
an order p = 2 (i.e., a square root), 


(7) 


Aan+1 ~ An 


a’ -N =0 

G, *@—)? h—0 

at + 2Zae, + e& —-N =0 
im — ae 
e, = ————- 
Dat e: 


or, ignoring the e, term in the denominator and _ setting 
Aant+1 = a F Ep, 


Aan+] =~ An + 


which is identical to Equation 6 for p = 2. 

To continue to Aitken’s fourth method, we can perform a 
second iteration of Equation 4. Since ag = 7 led to a value a, best 
represented as a fraction 50/7, it is useful to rewrite Equation 4 for 


a, = s/t: 
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t?-N — s2 


2s? 


An+1 = 4n + an ® 


Then, 


50 50 49 e 51 - 2500 
ee ee 
7 2 @ 2500 


ag 


50/7 + (50/7) © (-1/5000) 


TAAU4ZSD . .: 


which is Aitken’s fourth estimate. 

We can also develop an improved version of the Newton- 
Raphson method. This variation has been discovered many times 
and bears several names; we will adopt the more common and 
appropriate designation of Halley’s method, referring to Edmond 
Halley (1656?-1743) of comet fame [3,9,12—15]. 

In Equation 7 we can replace the (a,,; — a,) term in front of 
the second derivative f"(a,,) in the denominator with the approxi- 
mation found in Equation 5, which we obtained by ignoring this 
term completely. This is Halley’s formula: 


Lhhae)ot hae) 


a Ha len) fG,) 


For f(a) = aP — N = O, we arrive at a third-order formula for 


finding the pth root of N [6,16,17]: 


(p — l)ah + (p + 1)N 
(p + L)ah + (p — 1)N 


or, 
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+ | 


P 
ay + 


e (N - af) 
ane = Aq © (8) 


p +4 
ag + (1- ) oN —ap 
2p 


For the square root, we have p = 2 and 
2 + (3/4) # (N - a) 
ab + (14) © @N—-a) 


je) 


Aan+1 = an “ 


Returning to our example N = 51 and ao = 7, we obtain the 
first approximation in Aitken’s fifth estimate: 

ay = 7 (504/494) = 7.14141414... 

We can now rewrite Equation 8 for a!, = N/a?! as 


NP p+ ( a) 
i a. — ——= 


aP(p-l) 2p ap(p-!) 


NP p+ NP 
aP(p-1) 2p aP(p-l) 


After some algebra we arrive at the relation 


ae 
An+1 = an ® 


N?2-Paple-)) +(1—(p+ 1)/2p) e(N — N2-Pap(p-) 
N?2-Paple-)) + ((p + 1)/2p) ¢ (N - N?2-Pap(p-1)) 


' ane 
Aan+] = An 


(9) 


For p = 2, the extensive term multiplying aj, is the reciprocal 
of the corresponding term in Equation 8. This provides us with the 
second approximation by Aitken in the fifth estimate: 


al = (51/7) © (49%/50%) = 7.14142715700... 


Now Equation 4 gives a better approximation as the average 
of a, and N/a,. For a; = 7(504/49%), then N/a, = aj and we 


culminate in the fifth estimate, 
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ay = (a, + a})/2 = 7.14142842857 ... 


correct to eleven digits. 
Returning to Equation 8, we can transform the third-order 
relation to one that eases our later calculations: 


N — aP 
(p + 1)ah + (p — 1)N 


Ant] = an + day ® 


We now recognize this as second-order Equation 6 if we set N 
equal to aP in the denominator. 

Use of this algorithm requires nontrivial division when N is a 
multidigit number or p is larger than 2. Given this, we can write 
the relation in a form suitable for a, being a fraction s/t: 


fan) — g§P 
(p + 1)sP + (p — 1)t?N 


An+1 = an + da, $ 
or, for p = 2, 


t?N — s2 
3s? + t2N 


Aan+1 = ayn + dar, A 


Again, for a; = 50/7, 


50/7 + 2(50/7) a (60). 
= e@ — ~—_- — ~~ -— > 
F 3(50)2 + 49 © 50 
= 50/7 + (50/7) © (-2/9999) 

= 7.141428428557... 


which, since we are reducing 50/7 by one part in 49994, is 
apparently the extreme approximation used by Aitken in the sixth 
estimate to eleven digit accuracy. 
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The Chebyshev Correction 


We can add at this point a further generalization of the 
Newton-Raphson method discovered by P. L. Chebyshev 
(1821-1894). His paper on this subject earned a silver medal in 
a competition in 1840, but actually was not published until 1951 
[3,11]. His iterative algorithm is a third-order one, meaning the 
error falls off cubically with each iteration, or equivalently, that 
the number of correct digits is approximately tripled with each 
iteration. The relation is given by 


An+1] = An 


(an) (0) f"(a,) 
f(a,) \f(a,)/ — 2£"(an) 


Notice that this amounts to subtracting a term in each iteration 
from the Newton-Raphson algorithm. For f(a) = a? — N = 0, 


(Sony ; fap.) : (3 - =) ¢ p(p — Par 


fi(a,)/ — 2F (a) \ pak! 2pap-! 
a — gP\2 
_P 1 7 N — aP (10) 
dan \ aR 


For p = 2, this correction term is given by 
1 (~ - =f 

— °® 

8a, ai 

This method offers two advantages. It allows a preliminary 
estimate to be produced quickly using the Newton-Raphson 
method, and then improved afterwards with this correction if 
desired (merging the terms does not simplify the calculation). It 
also uses a divisor simpler than that used in Equation 8. 

For an initial estimate of aj = 7 for (51)!/*, we arrived at ay = 


50/7 by using the Newton-Raphson technique. We can therefore 
improve this estimate by finding 


Roots 87 


l 31 = (7)? 
a= 5077 - (5 )*( 7 


= 7.1413994 
compared to (51)! = 7.14142843 ... 


which is an improvement over 50/7 = 7.142857... but poorer 
than two iterations of the Newton-Raphson method, which earlier 
yielded 7.1414285.... As we will see for higher-order roots, 
however, two iterations of the Newton-Raphson method can in 
fact be prohibitive. For example, if a single-digit ap is used in the 
Newton-Raphson method to find a two-digit a, for the next 
iteration, then p doesn’t have to be very large to prevent the next 
iteration a) from being calculated, particularly since a? appears as 
a divisor. However, Chebyshev’s additional term uses the original 
single-digit ag in its calculation. 


A General Square Root Algorithm 


The selection of N = 51 is a convenient, though illustrative, 
one for Aitken and in general the approximations become unman- 
ageable very quickly. Incredibly, though, we can develop a method 
for mentally extracting square roots of more general multidigit 
numbers to greater and, with patience, even arbitrary precision. 

To begin, we rearrange the second-order Equation 6 for p = 2 
to express the amount b, to be added to a, in each iteration: 


WN - @& 


2a, 


Isaac Newton (1642-1727), incidentally, pointed out that hand 
calculations of square roots may be quickened if, after finding the 
square root the standard way to one-half the required number of 
digits, the remainder is divided by twice the existing calculated 
root to obtain the rest of the digits [4]. 

How good is this approximation? For our purposes, we can shift 
the decimal point (two places at a time) in the given number to 
scale the root to have two digits left of its decimal point. The 
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decimal point can then be moved in the root (one place at a time) 
when we are finished. For example, if N = 51, we can let N = 5100 
and find an approximate root 71.41, easily scaled back to 7.141. 
We do this because two-digit operations are most convenient in 
the general algorithm. 

We can further assume for the time being that a two-digit 
approximation ap to N!? jis achievable by hook or by crook. 
Remember that we can square two-digit numbers fairly easily. One 
convenient technique for finding ap is to begin with a round 


number. For N = 1867, 
1867 = 40 © 46 = 1840 

But 40 © 46 = (43)? — (3)2, so 
(1867)!/2 = 43 and 1867 — (43)? = 27-9 = 18 


Continuing now, 


18 
bo = = 9/43 = 0.20930... 
© 43 


and 
ag + bo = 43.2093... Actual: 43.2088... 
In general, the error in (ap + bo) is given by 


(N12 2. ag)? 


error = (ag + by) — N!? = 
Zag 

Obviously, for large ag the error is small, as —0.5 < (N!/? ~a,) <0.5. 
As it turns out, for a two-digit integer ag we find that 98% of the 
time the quantity (ap + bo) will be accurate to two decimal places, 
and 50% of the time to three decimal places, offering four or five 
digit accuracy in a square root with a minimum of effort. 

Since (ap + bo) will always be in excess, we can ask whether 
subtracting a number 0.000x before rounding, with x a digit, might 
improve the 50% average. In fact, it turns out that 0.0005 is the 
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optimum amount to subtract. Therefore, by simply truncating the 
result to three decimal places rather than finding the fourth and 
rounding, the result will be accurate to three decimal places 61% 
of the time. The third decimal place will be within one of the 
correct digit 85% of the time. 

The next iteration b,;, however, presents difficulties because 
a, = ap + bo is now a multidigit number; mentally calculating b, 
is generally prohibitive. We can in fact overcome this obstacle by 
developing an algorithm specifically for our purposes. 

Let ap and b,, be two-digit numbers defined by 


NI? = ay lag!" Tie... log... 


=aot+ }) 107%, 
k=0 
where the decimal place lies between ag and bp and the symbol 


“|b” as before indicates that digits in b are melded into the term 
to the left to leave a two-digit value of b. Now, 


k=0 


n—1 
ag & 2, 10h, 
k=0 


n—1 2 
IN - (a + y 10h, e 


10-"2p, = 


If we limit b, to two digits and approximate the denominator 
aS ap, we get 


n-1 2 

IN z= (1 4% 10h, |e 

10-2n-2p, = ii - li 
ao an 


where r,, is the remainder of the first term after dividing to two 
digits. Now, 
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IN _ (2 4 y 10%) | 


k = 0) 1 ieee 1 


eee SSS Se — 
ag ag 


but, 


n 2 
iH (a + » 10h, 
k=0 


n-1 2 
- (a +y 10h, + 107, 


k=0 


] 2 
i ot, 


-N-(a+¥ 


k=0 
1 


n-—- 
~ 2(10-2"-h,) « (a :e 107%, 

k=0 
— 10-442 


= 2610 ab. + 2010 "*) 


n-1 

— 2005") (a + y 107 hy 
k=¢ 

— 10a 


Therefore, 
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n—| 
ies, - 10>, )) 10-26, 10 yz 
k=0 
10-*4B...| - 
ag 
a ae 
ag 
Or, 
n—1 
100r, — b, ), 107*b, - 10-2"b2/2 
k=0 [n+] 
b,4.. = ——_—____—_ - 
ay ay 
2 
= genes IN = (@ol bo! br - ba) 


a 


This is perfectly general. The validity of the b,,; value is not 
affected by the approximations of ap in the denominators of earlier 
b’s. In fact, ay should be a rounded two-digit approximation to 
N!? so bo, which now may be positive or negative, has a magni- 
tude of 50 or less. 

Since we are concerned only with terms in the numerator that 
contribute to the division to two digits, we can shift other terms 
into the remainder term to be used in calculating b’s of higher 
order. Conversely, since each b multiplies the remainder of the 
previous b by 100, we need to extract these shifted terms in the 
b at which they become significant (i.e., when the power of ten 
preceding the term is 0). Therefore, for n even, 


%—1 
1OOR, — Ly baby - bi2/2 


k=0 
b+) = ——--—- + Rast 
ag 
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where R,, and R,,; represent the new remainders of b, and b,41, 
carrying with them the lower-order terms of previous b’s that have 
not become significant yet in the division to two digits. 


For n odd, 


yeaa 
100R4, = oc benibe 
k=0 
bas. = = 
ag 


Simply put, to find b,,, we subtract from 100R, pairwise 
multiplications of b starting from the outside (bob,,) and working 
inward to the center; if one b is left in the center, we subtract b?/2 
as well. We then divide by ag to two digits and use the remainder 
for calculating the next b. Despite the heavy mathematical 
notation and derivation, the concept is quite simple, and is best 
understood by working through an actual calculation. 

For example, to find the square root of 51, we first consider 
(5100)'", which provides a two-digit ao. 


70 © 72 = 4900 + 140 = 5040 

giving (71)? = 5040 + 1 = 5041 
Then, 

ay = 71 with a remainder of 59 
Root thus far: 71 


5900/2 
bo = eu = 4] Ro = 39 


71141 


We now begin the algorithm. 
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3900 — (41)*/2 3059.5 


b) = ———_- = —— = 3 R, = 6.5 
71 fol 
71141143 
650 — 41 © 43 a3 
»= = --(—) = ~(16R =23) = 216 
fi al 
R, = +23 
71141143 1-16 
y, = 2300 = (16) © 41 ~ (43)7/2 _ 2031.5 _ 
, 71 71 
R; = 43.5 
711411431-16128 
350 — 28 @ 41 — (-16) © 43 3890 
py, = £350 = 28 9 41 - (216) + 43 _ 3800, 
it 71 
Rg = 56 


711411431-16128154 
eee? 45 — (- 16)7/2 | 2054 — 


711411431-16128154128 


In the last step, we could have taken bs = 29, Rs = —5, but with 
the preponderance of positive b’s it seems likely that the numer- 
ator of bg will now be reduced to a number whose absolute value 
is less than 7100/2 = 3550 (giving be < 50), which is what we 
desire for two significant digit accuracy. In fact, the same philos- 
ophy was used in the two previous steps. If the absolute value of 
the numerator had turned out to be greater than 3550, giving 
a b = 50, we would have simply backed up and increased ot 
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decreased the previous b, adjusting the previous remainder accord- 
ingly. This is most likely to occur when ap is small or when cal- 
culating higher-order b’s involving several terms in the numerator. 

After all this buildup, it happens that I am mistaken and the 
numerator is still larger than 3550, but we’re close enough to call 
it good. 


6600 - 28 © 41-54 ° 43 + 16 ¢ 28 33Re 


711411431-16128154128150 

We could continue, but we’re already well beyond the eleven 
digits obtained by Aitken. Adjusting the decimal point back to 
where it belongs, we get 

(51)!? = 7.141428428542850... 

Actual = 7.141428428542850... 


An alternative in finding ay may be more convenient than what 
we did earlier: 


(70)? = 4900 


(5100 — 4900)/2 
70 


= Remainder = 30 


Therefore ag = 71 and 


=r 
bp ae Ro = 39 


and we are back where we were before. 
At any rate, we have succeeded in obtaining a sixteen-digit 
result with a reasonable number of two-digit multiplications and 
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divisions. In fact, all b’s in this algorithm have magnitudes less 
than 50 (except where we let it go in bg). Incidentally, did you 
notice that in by, for example, the terms (—28) © 41 + 16 ¢ 43 
can be reduced to (-12) © 41 + 32? This is useful in b; through 
bs. Two terms with multipliers of 28 also appear in the expression 
for be. 

Throughout this discussion I have been careful not to drop the 
divisor 2 into the denominator for a reason. If we alter the division 
by ag, we alter the remainder. Therefore, if the division is reduced 
by a factor m, we have to multiply the remainder by m. Also, it 
is incorrect to meld the b’s before the calculation ends. For 
example, the terms 431-16 cannot be converted to 42184 
when calculating b,’s, even though this is what is done in the 
final melding. 

It is also apparent that finding the square root of a four-digit 
number is no more difficult than for a two-digit number, since we 
expand the latter to four digits anyway. In fact, as in cross division, 
additional digits beyond the first four are simply taken pairwise, 
halved, and added to the numerator of the associated step. To 
demonstrate this, we will find the square root of 16460.89. Since 
(13)? = 169, we proceed from here: 


164 — (13)? = 5 


by =§ OE = UTR) = 17 Ro = 1 


131-17 


100 + 89/2 - (-17)2/2 
pe a so 


= 
i 
© 


We find in this particular example that (128.3)? = 16460.89 
exactly. Since b,; = 0 and R, = 0, we see that in no later step can 
the numerator be other than 0 as by = -17 will always be 
multiplied by a value b = 0. 

Therefore, we have a useful method that offers no significant 
increase in effort between finding the square root of a one-digit 
number and that of a multidigit number. The increase in complexity 
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with the number of root digits taken is quite low compared to that 
for the traditional square root method taught (at one time) in 
school. I encourage you to try extracting a sixteen-digit square 
root of 51 using the traditional method with paper and pencil; you 
will be astonished at the work involved. The integer-like character 
of this algorithm may also be useful for extracting square roots of 
extreme precision on small, integer-based computers. 

I emphasize at this point what may have slipped by without 
notice, namely that we are for all practical purposes outstripping 
the capabilities of pocket calculators and virtually all commercial 
computer programs. The main reason | ended our calculations at 
sixteen digits is the difficulty I had in obtaining a verification of 
the answer beyond it! This is a double-precision floating point 
value. It takes a specially designed computer program to exceed 
this, and these are not routinely available. In my view, we can 
excel over the available instruments in extracting square roots in 
two areas. First, we can generally find an answer to five digits in 
less time than a calculator can be retrieved and worked, much less 
a computer. Second, with time we can find an answer to extreme 
precision, beyond what a calculator or typical computer would 
have the capability of doing without special software. I think we 
can take some satisfaction in this. 

One expensive software package, a symbolic mathematics pro- 
gram, allowed me to arbitrarily define the precision of floating 
point values, and with it I obtained the square root of 51 to 
additional digits. The first 25 digits or so actually form an 
intriguing sequence of no significance that I am aware of: 
7.1414284285428499979993998 .... 

To my chagrin, but not unexpectedly, I found that the essence 
of this general method is given in other variations by Uspensky 
[18] and Lehmer [19], and perhaps others [20]. Lehmer derives the 
algorithm as a variant of cross division, run on the ten-digit 
mechanical computers of the day (1926). As in his cross division 
method, his algorithm contains internal checks on results and is 
therefore somewhat more involved than one designed for mental 
calculations. In addition, negative b’s are not considered, creating 
additional constraints on significant digits and requiring frequent 
multiplications by numbers of magnitude larger than 50. Also, 
since the divisor 2 appears in the denominator in determining 
quotients and remainders, the denominator generally consists of 


Roots 97 


one more digit than necessary. Uspensky traces the algorithm to 
Fourier Division using single digits and again does not consider 
negative b’s. 

To implement this algorithm using single-digit values of b, 
though, is more tedious and, in the end, problematic due to poorer 
resolution, as “backing up” can affect twice as many b’s. Uspensky 
discusses a convenient correction term when multiple b’s are 
altered. Three-digit values lead to frequent reversals to signifi- 
cantly adjust previous quotients. Two-digit values provide a nice 
compromise conducive to mental calculation. 

In conclusion, I mention that for extracting the square root of 
a fraction, it is usually better to obtain a denominator that is a 
perfect square: 


(2/7)! = (14/49)! = (14)!?/7 


Beyond its immediate use, the ability to readily extract square 
roots offers an enormous advantage in developing other algo- 
rithms, including factoring methods that require an integral square 
to be identified, approximations of roots of higher order, very 
precise trigonometric approximations, and so forth. It is a very 
handy tool to master. 


The Reciprocal Square Root 


The reciprocal square root, or N~!””, can be calculated as a 
special case of the Newton-Raphson method applied to the 
equation N — l/aP = 0 [6,21]. Using Equation 5, we arrive at 


an 
Sst -~a*_— * (1 — Nab) 


which we can write for p = 2 as 
aye) =, + Say(l — Ned) 


a relation surprisingly free of divisions. 
For example, let us find the value (51)~!/? = 0.1400280084028 . . . 
(a number, incidentally, as fortuitous for us as 51 was for Aitken). 


We can take ag = 0.14. Then 
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ay — ag = .9(.14)[1 — 51(.0196)] 
= .07(.0004) 
= .000028 
ay = .140028 


Furthermore, multiplying the above result by N provides an 
approximate value for N!”, obtained without significant division: 


51(.140028) = 14.0028/2 + .140028 
= 7.141428 
Applying Halley’s method to the equation N — 1/a? = 0 results in 


F 2a,(1 — Nab) 
2p — (p + 1) © (1 — Nab) 


An+1 = An 


Again, for p = 2, N = 51, and ag = 0.14, 
a; = .14 + 1.12(10~)/3.9988 
= 1400280084023... 


This is a stunning result, but the division is complicated. 
Perhaps the best approach is to factor out a 4 in the denominator 
and approximate 1/(1 — .0003) as (1 + .0003). Then a, = .14 + 
2.8(10-°) © (1 + .0003) = .1400280084. 

Another alternative, generally more advantageous when ag is 
less accurate to the true value N~!?, is to derive Chebyshev’s 
correction term for f(a) = N — l1/a? = 0: 


(Sey my) JON - 2°)’ - pie + dae 
ae fee”! | LEE eee 
ff (ae) 2f (a,) p-a2(p+ 1) 2pa-(P* 1) 


Eee i! 
op 


ap(1 — Niap)- 
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For p = 2, N = 51 anda, = .14, we can subtract this term from 
our earlier result from the Newton-Raphson method: 


a, = .140028 + (.14)(.0004)? 


2° 4 


.1400280084 


accurate again to the last digit. 


Cube Roots 


Although the concept used in our last most general square root 
algorithm does not apply to extracting higher-order roots, we may 
revert to Equation 6 as a workable second-order approximation. 
We can rewrite this as 


(p aaa 1) an, + N/ak! 
Pp 


In other words, a,,; is the weighted average of a, and N/aP!. 
For this relation, the error e; in a, is less than (p — 1) eo/p. In 
general, the error in the nth iteration, for e; the error in aj, is less 


than [22-24]: 


(p — 1) ey 
dans 


forn> 1 


Thus, after a; is found, the error in each step diminishes as the 
square of the previous error, defined as a second-order relation and 
approximately doubling the number of correct decimal places in 
each iteration. 

As an example, consider (119)'3 = 4.9186847..., which we 
initially approximate as a) = 5. Then 


gis + 119/(5)? 


= 4.92 
3 


aj 
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Choosing ap = 4.9 provides an a; = 4.91875 ..., correct to four 
digits, but the division is more difficult. 

Returning to the first weighted average a, = 4.92, we can 
increase our accuracy by subtracting Chebyshev’s term, given in 
Equation 10, for p = 3: 


1 119 — (5)\? 
4.92 — . py 
905 (5)2 


= 4.91872 


This result is more accurate than the Newton-Raphson method 
with ag = 4.9, and was obtained with much simpler division. 

An alternate approach for cube roots involves, to gratefully 
adopt Aitken’s terminology, the process of “thirding,” a method 
corresponding to the quartering performed earlier for the square 
root [1]. Equation 8 provides a general method for the pth root, 
and for p = 3 it becomes 


Bee) © (IN - a2) 
a+ (mey © (N’- ag) 


Aan+1 = a, ° 
Again choosing ap = 5 for (119)!/3, we find (N — ag) = -6 and 
ya 
a, =5 © — = 49186992... 
123 


which is a more accurate result than before, although this partic- 
ular example involves a three-digit divisor. 

The complementary process of “sixthing” operated on ap = 
119/(5)* is another technique. Equation 9, the basis for the 
corresponding square root approximation, proves intractable for 
p > 2, so we turn again to Equation 8 with a,,, = N/aP3}: 
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a (- ak + ((p + 1)/2p) © (N - af) " 
P+ (1 - (p + 1)/2p) © (N — af) 


If we define y = 4(N + ab), 


i—* — 
ak! \y + (1/2p) © (N - a?) 


N — 1 E BL 
ap! 2p 8p? 


—] 
x yP O(N — ab)? +.. ) + (se - 7 y?-2(N — aP) 


# = a yP 3(N — ab)? +.. ) 
8p? 


by the Binomial Theorem (discovered by Newton). 

Now since y is the average of N and af and (N — af) is assumed 
small in comparison to y, the terms after and including (p — 1) x 
(o— 2) 9? - aP)?/8? can be eliminated with minor error, leaving 


ap-! = i 
: y??(N — af) 


Se. oe (11) 
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For p = 3, we have the sixthing process: 


Na + (1/6) © (N — a?) 
ae = 8 


ab ap + (iG) © (N - a) 


Why would anyone use sixthing when any (N — a2) divisible 
by 6 would also be divisible by 3? The answer is that sometimes 
the numerator and/or the denominator is simpler in sixthing, and 
also for very small error in a,.;, the weighted average of the two 
techniques provides even greater accuracy. 

Indeed, N = 119 was chosen to demonstrate the first advantage. 
For a = 119/(5)?, 


ay = (119/25) © (124/120) = 4.76 + 4.76(1/30) 
= 4.9186667... 
Aitken remarked in 1954 that he was not aware that sixthing 
had ever received notice before, although he gave no basis for it. 
We can expect marginally better accuracy (perhaps to the fifth 


decimal place) if we find the weighted average of a, and a}: 


' 
day * & 


= 4.9186 + 10-4(1/3) © [2(.99) + .67] 
= 4.918688 ... 


an approximation accurate to six digits. 
Again, an alternate cube root formula [25] can be derived from 


e — N!/3)3 = aa o Mz 3a2N13 A 3a,N23 
Now an approximation a,,; satisfies 
at.1 — Aanaaei + (a, — N)/3a, = 0 


giving, for 0 < a3 < 4N, 


4N — ad\l2 
an+1 = 4 a + ( ) 
Ga. 
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The square root term to several places is manageable (using 
methods given earlier) and the formula is of third order. 


For N = 119 and ap = 5, 


476 — 125\12 
a = {5 + ( Tes | 


= ¥ [5 + (351/15)"] 

= ¥ [5 + (23.4)!?] 

= (5 + 4.837355...) 
= 4.918677... 


again accurate to six digits (and better here than thirding or 
sixthing alone). In general, this method offers the advantage of a 
smaller divisor, 3a,, compared to roughly a2 in thirding or sixthing 
(unless simplification is possible, as for N = 119). The obvious 
disadvantage lies in evaluating the square root. However, the 
insensitivity of the square root algorithm to the number of digits, 
as opposed to division by a multidigit number, may quickly propel 
this method into the forefront when a, contains two or more digits. 

If we really don’t mind taking a square root (as by using the 
general square root algorithm given earlier) we can achieve even 
greater accuracy in a cube root without raising a,’s to powers and 
with division limited to that within the square root extraction 
process [26]. To do this, we perform the Newton-Raphson proce- 
dure not on the equation a’ — N = 0, but on the equation 
aj? _ N!/2 = 0 or a3/4 — N/4 = 0. We then produce 


an+1 = (1/3) © [2(N/a,)"? + a] (12) 
and 
ans) = (1/3) © [4(Na,)!* — a,] (13) 


both more accurate than the original Newton-Raphson iteration. 
For N = 119 and ap = 5, Equations 12 and 13 yield a, = 
4.91901624 ... and a, = 4.91851830. .. respectively, compared to 
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the actual value of 4.91868473 .... For ap = 4.9, we obtain a; = 
4.91870254...and a; = 4.91867584 .... However, notice that if 
we can manage to multiply N and a, in Equation 13 for a highly 
accurate a, or at least divide N by a, in Equation 12, we are on 
the way to an extremely accurate (if time-consuming) cube root. 
This assumes that we are comfortable taking square roots to 
several places (and I would rather extract a square root than 
perform a long division). 

For example, we can find 119(4.92) = 119 ¢ (5 — .08) = 585.48, 


and from Equation 13, 
a, = (1/3) © [4(585.48)!/4 — 4.92] 
= 4.91868469... 


accurate to eight digits. 

While not applicable to N = 119, we can sometimes recognize 
that by multiplying or dividing N by the cube of a small number, 
or the ratio of cubes of small numbers, we can arrive at a number 
close to a cube. For example, for N = 91, a number distant from 
the nearest cubes (64 and 125), we can multiply N by 8 to arrive 
at 728, very close indeed to 729, the cube of 9. Therefore, N!? 
is approximately 9/2 = 4.5, where 2 is the cube root of 8. Better 
yet, we can take a weighted average: 


728/(9)? = 8.98765... 
2N!3 = (2 ¢ 9 + 8.98765... )/3 
N!3 = 3 + 8.98765/6 
= 4.4979424 ... compared to 4.4979408 ... 

Failing this, we can often divide a number by 2 or 5 to produce 
a new number close to a cube. We then find the approximate cube 
root of this new number, and adjust the result by multiplying by 
(2)'8 = 1.25992... = (1 + .26) or by (5)! = 1.70998... = (1 + 


71). The factor (3)! is less convenient, and is perhaps best 
approximated by the fraction 10.1/7. 
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Higher-Order Roots 


Since fourth roots can be extracted as two successive square 
roots, we continue on to fifth roots, say to (48)! = 2.16894354.... 
In general, a second iteration of any technique is prohibitive, so 
we can explore algorithms that increase the complexity of each 
iteration. For an initial estimate of ay = 2, a weighted average of 
ay and 48/(2)* gives a; = 2.2. Subtracting Chebyshev’s term results 
in a; = 2.16. Three-fifthing, from Equation 8, produces from ag the 
value 2(41.6/38.4) = 2(13/12) = 2.1667.... The complementary 
process of tenthing, from Equation 11, produces from ag the value 
(48/16) © (33.6/46.4) = 3(21/29) = 2.1724.... Their weighted 
average yields 2.1678 .... 

Another method for higher roots like this is to write the 
equation N — aP = 0 as 


a a oy 
he 


Now, for —a, << (a —a,) << a,, we can use the relation [27,28], 


fi leg 
— ee + —__— 
n q ‘an q 


to convert the equation into one of a degree that simplifies the 
process of solving it. Substituting, we get 


This can be iterated: 
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If we stretch this method to the limit, we can let the new degree 


q be 1; then, for p = 5, N = 48 and ap = 2, 
a, = (1/5) © 48/(2)* + (4/5) © 2 
= 


This is the same answer we obtained from the weighted average 
earlier, and indeed this relation always reduces to a weighted 
average if q = 1. (Incidentally, cascading the equation down to 
q = 1 by steps of 1 results in the same formula as above for aj). 
Therefore, we can expect to do better if we reduce the order p by 
less. For example, we can let q = 4 and solve for a; by successive 
square roots: 


af = (4/5) © (48/2) + (1/5) © (2)4 = 22.4 
at = 4.7328... 
ay Sa 2.1759 eae 


a procedure still not as accurate as two-fifthing or tenthing, much 
less their weighted average. However, since we only divide here 
by ap, not by powers of ap as in the previous algorithms, we can 
stretch our capabilities of multiplication and actually use a; = 2.2 
to calculate an iteration: 


ad = (4/5) © (48/2.2) + (1/5) © (2.2)4 


Now in this case, we may recall that the pth power of 11 is the 
single-digit melded (p + 1)th row of Pascal’s Triangle. The fifth 
row is 1,4,6,4,1, so (2.2)* = 16(1.4641) = 23.4256. Alternatively, 
since (22)* = 484, then we find (484)? = 500 e 468 + (16)*. Or 
again, we can use the rule of 25 on (48.4)? = (25 — 1.6) 100 + (1.6)? 


and remove the decimal point. In any event, we soon arrive at 
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aj = 22.139665... 
as = 4.7052806... 
a, = 2.169166... 
our closest approximation yet to 2.1689435.... 


Yet another approach for fifth roots was published by DeLagney 
in 1692 [3,14]: 


54 byt ((2 b =). =i a a 
a ={{—+—]} -— — 
a 4 4 2 


For a = 2 and b = 16 (not, as you might imagine, an ideal 
circumstance), we find 


(48)! = [(1.6 + 4)? - 1)? + 1 
= (13664319 ...)!? + 1 
= 2.1689448 ... 


which is extremely good, considering that b is not small compared 
tora’. 

Further yet, it is not wholly unreasonable at this point (or is it?) 
to perform this approximation with a = 2.2. As before, we arrive 
at a* = 23.4256 and therefore a*/4 = 5.8564. Now 2(23.4256) = 
46.8512 and a = 46.8512 + 4.68512 = 51.53632, giving us the 
result b = —3.53632. Then, 


(48)"5 = [(-3.53632/11 + 5.8564)! — 1.21]'? + 1.1 
= (5.5349164...— 1.21)! + 1.1 
= 2.16894364 ... compared to 2.16894354... 
In addition to this method’s higher accuracy over others, it 


offers very simple divisors in its terms, an important benefit for 
mental use. Halley generally preferred these approximations 
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involving square roots to rational approximations (fractions whose 
numerators and denominators are polynomials) because “. . . mani- 
fold Experience has taught me” that square root extraction is 
easier than division by a large number [14]. 

An alternative is to use an approximation, again by Halley, that 
contains only one square root extraction [14]: 


p-2 
a + 


i [a’ + 2ba??(1 — 1/p)]!” 
— p — 


(aP + b)1/P = 


which falls directly out of a new iterative relation for finding real 
roots of a real polynomial equation: 


Mag) ag)? = Blog) an 
f"(a,) f"(a,) 


An+1 = Ap 


Now for p = 5, 
(a? + b)!? = 3a/4 + (1/4) © (a? + 8b/5a?)'” 
Repeating the above example of a = 2, b = 16, we find 
(48) = 3 © 2/4 + (1/4) © (4 + 16/5)!” 
= 1.5 + (2166928 .. . )/4 
= 2.1708... 


This is not too bad considering the reduction in effort over the 
formula involving nested square roots; on the other hand, the case 
a = 2.2 is prohibitive here. 

You don’t have to do many of these approximations, particularly 
for numbers far from a perfect power, before it becomes apparent 
that, in general, the way to extract higher-order roots is through 
logarithms. The conversion to and from logarithmic form is 
not as hard as it may seem, and we will address this subject in 
the next chapter. 

How far can we go, then, in mentally extracting roots? At least 
for square roots, we can with our general method achieve a result 
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whose accuracy depends only on our patience and care. The 
mathematician John Wallis (1616-1703) wrote in 1669 that “In 
a dark night, in bed, without pen, ink or paper or anything 
equivalent,” he mentally extracted the square root of (in 
his notation) 


3,00000,00000,00000,00000,00000,00000,00000,00000 


and found it to be 1,73205,08075,68877,29353 to 21 digits. Two 
months later he wrote that he did the same for the 53-digit 
number 


24681357910121411131516182017192122242628302325272931 


and found the result to be 157103016871482805817152171 to 
27 digits [2]. 


I think we can take this as something of a benchmark. 
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Logarithms and 
Their Inverses 


It can be said without exaggeration that our modern 
cities with their bridges, railroads, factories, power 
installations, skyscrapers, in short, all our highly 
artificial environment, built by mankind during the 
last century or two, and which implies an enormous 
amount of numerical computations, would be 
impossible without the tables of logarithms. 

E.G. Kogbetliantz (1968) [1] 


Vine principles of logarithms were discovered first by Archi- 
medes, then rediscovered and used by Iranian astronomers in the 
thirteenth century, and finally discovered again and published by 
John Napier (1550-1617) in 1614. In addition to spontaneously 
arising as functions in the solution of mathematical and physical 
differential equations, they provide us with tools for vastly sim- 
plifying arithmetic operations such as multiplication, division, 
powers, and roots. In fact, Napier’s original intention for these 
functions was to ease calculations of trigonometric formulas, 
demonstrating their range of applications. Specifically, when nau- 
tical tables were unavailable or simply not used, the trigonometric 
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calculations needed for estimating position and time were histor- 
ically done using logarithms. Slide rules were strictly based on 
logarithms, and many computer routines rely on them as well. 

The logarithm of a number N to a base B is defined to be the 
value x such that 


BX = N 
Therefore, 
x = logpn N 


where the subscript B denotes the base of the logarithm. “Com- 
mon” logarithms have a base of 10, and the notation log N with 
no subscript is usually considered, as it is here, to have this 
decimal base. The common logarithm x = log N has an inverse, 
then, given by N = 10%. 

A base not tied to the base of our number system, but providing 
a convenient, or “natural,” base for logarithmic solutions to 
differential equations, is denoted by the letter e. The value of 
the transcendental number e in our base 10 arithmetic is 
2.718281828 ..., and the notation for the logarithm to base e, 
the natural logarithm, is generally given by In N. The inverse of 
the natural logarithm x = In N is given by the exponential 
function N = e*. 

Logarithms are useful in arithmetic calculations because they 
have the following properties: 


log (ab) = log a + log b 
log (a/b) = log a — log b 
log (a) = b log a 
log (a!) = (log a)/b 
These are true regardless of the base of the logarithm. 
In short, logarithms reduce multiplications and divisions to 


additions and subtractions, and convert powers and roots to 
multiplications and divisions. For our purposes, their latter use is 
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of the most benefit. To find, say, the seventh root of 13781, we 
could find log 13781, divide by 7, and find the answer as 10 to 
that power. The accuracy of the answer is limited only by the 
number of digits retained in the calculation. Also, it may be easier 
to find either log N or In N, depending on N, but it is almost 
always easier to raise 10 to a power than to raise e to a power. This 
ultimately requires us to work with common logarithms when 
inverses are required. 

To find a means of converting from common logarithms to 
natural logarithms (and back), we begin with relations that, by the 
definition of a function and its inverse, are given by 


10 = eln 10 
e = 10Q!°e e 

Substituting the second into the first relation, we have 

10 = (10!ee e)in LD — 10!°s e In 10 
Therefore, log e In 10 = 1. For a number N given by 

10!°2 N — x= eln N 
we again substitute the relation e = 10!°2 ¢ to get 

1Q!es Wo 10!°s e In N 
or, 

log N = log e InN 

Therefore, we arrive at log N by multiplying In N by the 
constant log e, and we use the reciprocal of this constant to 
convert from log N to In N. From formulas given later, we can 
find log e = 0.4342945 ... (or about .4343) and (1/log e) = In 10 = 
#302565 ..., so 

log N = .4343 In N 


In N = 2.303 log N 
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Again, we often are not interested in finding inverse logarithms; 
natural logarithms appear very often in formulas requiring 
calculation. 

In any event, the conversion of a number into logarithmic form 
and back is a useful and entertaining knack, and the tools for 
mentally doing so are the subject of this chapter. 


General Logarithmic Approximations 


At first glance, it may appear to be an enormous task to 
calculate logarithms, as it was for Napier, Briggs, and others. We 
see first that log 1 = 0 (a definition, actually) since (10)° = 1. Also, 
log 10 = 1, log 100 = 2, and so on, but intermediate values offer 
little hope. We can take, for example, (10)!/? = 3.1622777..., 
but this just gives us log (3.1622777 ...) = 0.5, something akin 
to throwing darts at the problem. 

The general procedure, and one deduced by Henry Briggs 
(1556-1630) in devising his tables of logarithms published in 1617 
and 1624, is to use the properties of logarithms given earlier to 
reduce the argument of the logarithm to a value very close to 1. 
Then, truncating the following well-known power series provides 
the natural logarithm to the degree of approximation desired: 

Zz x? x4 
In (1 +x) =x-—+—-—t+... for-l<x<1 (14) 
- 3s & 


It is the cleverness we exhibit in reducing the logarithm to a value 
near 1 that makes calculating logarithms a most creative and 
engaging activity, and lightning calculators have had a deep 
appreciation of this [2]. 

A very interesting description of Napier’s and Briggs’ develop- 
ment of logarithmic tables, as well as the efforts of others later, 
is found in Goldstine [3]. We will consider these and other avail- 
able techniques best applied to mental calculation. 

Obviously, the optimum procedure is to leave x = 0, which 
amounts to factoring the argument into powers of ten and low 
primes whose logarithms are memorized. If that is not possible, we 
can factor a number very close to the argument and extract a value 
of x for the power series correction. To illustrate, 


i 
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log 1200 = log [3 © (2)? © 100] 

=log3 +2 ® log2 +2 

seque2 + 2(.30903) + 2 

= 3.07918 
log 1211 = log [1200 © (1 + 11/1200)] 

= log 1200 + (.4343...) © [11/1200 

=/7A(iijigeey? +...] 


Taking only the first term of the series expansion and letting 
.4343 ... become simply .43, we have 


log 1211 = 3.07918 + .43(11/1200) 
= 3.08312 Actual: 3.08314 


If we had preferred, we could have divided by 2.302585 = 2.3 
instead of multiplying by .4343. Dividing by 2.3 gives log 1200 = 
39083117. 


Alternatively, we could have taken 
log 1211 = log 1210 + log (1 + 1/1210) 


since 1210 is easily factored into low primes and a power of ten 
and the value of x is smaller. 
A variation on this was used by the mental calculator George 


Parker Bidder (1806—1878) [2]: 
log (1 + x) = 10™ © x @ log (1 + 10°™) for 10O™< x < 10™"! 


The value of m is simply that which makes the term 10™x lie 
between 1 and 10. We need to know log (1 + 10°™™) as well: 
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log 1.01 = .00432... 
log 1.001 = .000434... 
log 1.0001 = .0000434... 
etc. 
Obviously, as m increases, the digits of log (1 + 10-™) approach 
those of log e = 0.4343.... 
To repeat the previous example, 
log 1211 = log 1200 + log (1 + 11/1200) 
= 3.07918 + (11/1.2) log 1.001 
= 3.07918 + .00398 
= 3.08316 


Where does this approximation come from? We know from 
Equation 14 that 


In (1 + x) = x 
In (1 + 10°") 2 10™ 


In (1 + x) _ log (1 + x) 


= (x 
lal + 10e) Wei + 10) 


log. {1 +30) TO "eo ¢ log.(1 + 10™) 


Now, since m is chosen so that x and 10-™™ are about equal, the 
errors in truncating their series are of the same order; therefore, 
this approximation is more accurate than that obtained by simply 
truncating the series for log (1 + x). 

It becomes clear that we can vastly improve our performance 
if we memorize the logarithms of a range of low-order primes, say, 
up to 11. 
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log 2 = .30103... 
lag 3 = 47712... 
log 5 = log (10/2) = 1 — log 2 
log 7 = .84510... 
log 11 = 1.04139... 


Actually, it’s an intriguing pursuit to estimate these as well, and 
very useful if they’re not memorized [4,5]. For example, 


(2)!0 = 1024 = 1000 
SO, 

10 log 2 = 3, or log 2 = 3 
Here we could really improve the estimate: 

10 log 2 = 3 + .4343(24/1000) 

log 2 = .30104... 

and if we dare to include the next term in the series expansion, 
~.4343(.024)*/2 = -.0001251..., and divide by the factor 
10 again, 

log 2 = .3010298... Actual: .301029996. .. 
illustrating that we don’t have to have a large number to arrive 
at a small value of x if we use our imagination. 

Other relations we can come up with (and improve with 
correction terms) include: 


34 = 81 =23 @ 10 
34 © 5? = 10125 = 104 
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7? = 49 = 10°/2 
3 0 2 = 1029 = 10° 
24 @ 37 e 7 = 1008 = 10° 
117 =i 2S 38 fp 
7 ¢@ 11 © 13 = 1001, or log 13 =~ 3 — log 11 — log 7 
27 © 37 = 999, or log 37 = 3 — 3 log 3 


The last two products we recall from the Chapter 2 discussion of 
divisibility (notice the result if we multiply them together!). In 
fact, one reason for the amount of time spent earlier on factoring 
numbers is the application here in calculating logarithms. 

We can create many of these relations, some very complicated 
but precise [6]. One type of relation that is generally useful for 
finding the logarithm of a prime p is 


p’=(p+i1)¢(p-1)+1 
=p + 17 * ap 7) 


All prime factors of (p + 1) © (p — 1) must be less than p. 
Mermin [4] gives the following examples (there are obviously 
many more): 


(37 e 11)? = 100 © 98 
74 = 50 © 48 = (100/2) © 3 @ 24 
217 = 22 © 20 = 2? © 11 ¢ 10 


We can even use this in our example of finding log 1211. In 
this case, the error in dropping 1 is extremely small. 
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201) = 122 1210 + 1 = 1212 ¢ 1210 


= 12 ¢ 101 e 121 @ 10 


=2?e3¢10¢ 11’ 101 


log 1211 = log 2 + log 11 + % [log 3 + log 10 + log 100 
+ .4343(.01)] 
= .30103 + 1.04139 + % (.47712 + 1 + 2 + .00434) 
= 3.08315 
Rounding errors produced the error in the last digit. 


Neighboring Value Relations 


Another form of relation for the natural logarithm In N (which, 
as earlier, is easily converted to the common logarithm log N) uses 
a variable transformation [1,7]: 


» fe! 1 +u 
N + 1 au 


u 


Then, from Equation 14, 
In N = 2(u + w/3 + w/5 +...) (15) 


where we have now eliminated the even terms in the expansion. 
However, this is mostly useful for N near 1, as the series converges 
very slowly otherwise. 

A much better approach in general is to take Equation 15 
and replace N with the expression (N + 1)/N. Then u becomes 
1/(2N + 1) and Equation 15 reduces to the very nice result: 
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2 
In (N+ 1) =InN + 
2N + 1 


«1+ + | (16) 
3(2N + 1) 5(2N + 1)4 


This converges very rapidly for larger values of N. 
Notice that if we take one term of the earlier series expansion 


(Equation 14), we find 
In (N + 1) = In [N(I + I/N)] 
= InN + In (1 + 1/N) 
~ In N + 1/N 


The new series (Equation 16) produces a more accurate value 
when the first term is taken, and converges more rapidly: 


2 
In (N + 1) =InN + 
ZN + | 


Therefore, setting the conversion factor .4343 ... to simply .43 for 
this small second term, 


log 1211 = log 1210 + .43(2/2423) 
but, 

log 1210 = 2 log 11 + log 10 = 3.08279... 
sO, 

log 1211 = 3.08279 + .00036 


= 3.08315 
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Obviously, more digits are needed here to show the greater 
accuracy, a pleasant situation after all. Of course, we are consider- 
ing here precision to the fifth decimal place; in reality, we can 
reduce this precision if we want to simplify the calculations. 

We can also replace N in Equation 15 by (N — 1)/N, giving us 
a relation for In (N — 1): 


la (N — 1) "In N - 
mw -— | 


ee | 
2m t) 5@QN - 1) 


This is not as useful for finding log 1211, as 1212 does not factor 
quite as nicely as 1210 does. 
In general, if a is a positive or negative number, 


2a 


la GN +a) = In N - 
2N +a 


] 
x|1 ++»... 
3(2N + a) 5(2N + a)4 


Intermediate Value Relations 


There is also a nice formula due to Halley for finding the 
logarithm of a number N equidistant between two numbers a and 
b whose logarithms are easier to calculate [3]. For N = A(a + b), 


log (ab)'" = log aia” 

/(a + b) 

4ab 
(a + b)? 


= /, log 


(b - a 


= ¥ lo 1 - 
2 10g i + i 
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or, from Equation 14, 
log (ab)! — log N = 4% (.4343...) 


y| ra ,, (b - a) 


Es oy P a tas by ie 


For (b — a) small compared to (a + b), this relation converges 
very fast due to the squaring of the denominator. Retaining only 
the first term of the expansion, we arrive at 


(17) 


(b — =, 


| N= 411 log b + .4343 
og N =X |log a + log b + 4343 ———— 


Now for b — a = 2, we can derive an amazing formula in the 
following manner: 


4ab , 4a(a + 2) 
(a+byY fae G + op 


f 2a* + 4a 


2 


j es + 44 + 1) + |. 
(2a2 + 4a + 1) -— 1 


. E + 1/(2a’ + 4a + a 
1 — W@a* + 4a + 1) 
If we make the transformation 
y’ = 2a’ + 4a +1 


(where the square of y is denoted to indicate the order of 
magnitude of the result), then 
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(2b) tog Lt 
ee Ale 


— log 


or, from Equation 15, 


(ab)!/2 


| l 1 
— log -(4$343...) 0/544... 
by 


3y° 5y!° 
It is difficult, I think, to appreciate the rate of convergence of this 
expression without an example. Halley provides an example of 
finding log 23: 
y= 1057 


(22 © 24)! 
SS — 


= (.4343...) - 


—— + + 
1057 3542796579 


Even for small numbers, then, the second and later terms are 
extremely negligible. Continuing, 


log 23 = % (log 22 + log 24) + .43/1057 
log 22 = log 2 + log 11 = 1.34242... 
log 24 = 3 log 2 + log 3 = 1.38021... 
average = 1.36132... 
log 23 = 1.36132 + .00041 


= 1.36173 


This result is accurate to the last digit, which is all the accuracy 
we were retaining in the calculation. The formula is obviously 
capable of extreme accuracy, however, and for those who anti- 
cipate pursuing it, a convenient multiplier to replace .43 in 


126 Dead Reckoning: Calculating Without Instruments 


approximating .4342945...is given by the fraction 43/99 = 
434343... (the value 1/2.3 = .43478...). The best alternative 
we have found so far in calculating the logarithm of a relatively 
low prime number such as this is: 


(23)? = 22 © 24+ 1=22 © 24 
log 23 = 4 (log 22 + log 24) 


which is seen as simply the first term of the last result. 

Halley’s method does require division by y? and multiplication 
by .4343..., although we only kept two digits in these opera- 
tions. However, we can eliminate even these nuisances, given that 
we know the logarithms of the neighboring numbers, at the cost 
of some of this extreme accuracy [3]. This relation is due to the 
ubiquitous Newton and for practical purposes the reduction in 
accuracy is not important. 


Newton defined for a = N — x and b=N + x: 


d = % (log b — log a) 


Then, 
W + x 
2d = log 
— x 
1 + x/N 
= log 
1 —- XN 
and from Equation 15, 
xX oo x? 
2d = — + — 
N 38 SNP 


Also, the difference (log N — log a) can be denoted as g: 


N 
e- Ge = — log (1 — x/N) 
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From Equation 14, 


x2 x? 


0256 = 
» 2 OW 


Therefore, 

x x? x? 
es N 2 3h 
ta x x? x 

— a ee See 

N 3N? 5N>? 


Dividing through, we arrive at 


%dx (1/12) dx? 
+ — + ————_ + 


=d = 


OT, 


q dx 
rr S| a 
, N 


Since g is the amount to be added to log a and d = % (log b 
— log a), we can rewrite this as: 


/Ad 
log N = %4 (log a + log b) + —— (18) 


As d implicitly contains the constant .4343...in finding the 
difference between common logarithms, we are rid of this multi- 
plication. In addition, we divide by a number N, not by the much 
larger value y’ = 2a” + 4a + 1 of Halley’s method. It is also easier 
than Equation 17 for x > 1. 

Returning to our previous example for log 23 with x = 1: 
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| Yd o | 
log 23 = 4 (log 22 + log 24) + 


log 22 = log 2 + log 11 = 1.34242... 
log 24 = 3 log 2 + log 3 = 1.38021... 
average = 1.36132... 
d = .018895... 
log 23 = 1.36132 + .00945/23 
= 1.36173 
The result is identical to Halley’s method to the accuracy 
retained throughout the calculation. Notice that we can simply 
subtract log 22 from the average to find d, since 4 (log a + 
log b) — log a = % (log b — log a). Now if we actually knew 
log 22 and log 23 to extreme accuracy (and multiplied by 
4342945 ...in Halley’s method), we would find these results: 
Halley’s Method: log 23 = 1.36172783590... 
Newton’s Method: log 23 = 1.36172770649... 
Actual Value: log 23 = 1.36172783602 ... 
The impressive accuracy of Newton’s method for x = 1 leads us 
to try log 23 for x = 2. While perhaps necessary in general for 
some numbers N, the case x = 2 is actually a little easier for 


log 23 as well if log 7 is known. 


| Yd © 2 
log 23 = 4 (lop 21 + loge25) + 


log 21 = log 3 + log 7 = 1.322219... 
log 25 = 2 log 5 = 1.397940... 


Logarithms and Their Inverses 129 
average = 1.360080... 
d = 03786... 
log 23 = 1.360080 + .03786/23 
= 1.36173 


Again, this is accurate to the number of digits retained. More 
exact values would give log 23 = 1.361725754.... 

In short, Newton’s method gives highly accurate results for 
reasonable values of x without the more difficult multiplication 
and division of Halley’s method. 

This leads us naturally to consider another more general inter- 
polation technique. For a known set of values f(x), f(x;),..., 
f(x,,) located near or around the unknown value f(x), where xo, 
X1,-..-, X, are not necessarily equally spaced, the most straight- 
forward technique for mental calculation of f(x) is through 
Lagrange’s formula [3]: 


I] (x — x,,) 


nN ok#m 


f(x) = ©” 
M0 Lia 
k#m 


f(x;) 


where again the Greek letter sigma here represents the sum of 
terms for k = 0 ton, and the Greek letter pi represents the product 
of all terms except for k = m. 

Let us consider the calculation of log 13 = 1.1139434.... 
Since log 11, log 12, log 14 and log 15 may be easily generated 
from factors (assuming the logarithms of primes < 11 are memo- 
rized), we would expect a better result than Newton's formula 
(Equation 18), which involves just log 12 and log 14. 

For x, of this spacing, Lagrange’s formula reduces to: 


* 2 u | 
log 13 = — — log 11 + — log 12 + — log 14 — — log 15 
og Pe | 3% he , me 


log 11 = 1.04139... 
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log 12 = 2 log 2 + log 3 = 1.07918... 
log 14 = log 2 + log 7 = 1.14613... 
log 15 = log 3 + (1 — log 2) = 1.17609... 
log 13 = (1/6) @ [- (1.04139...) + 4(1.07918...) 
+ 4(1.14613...) — (1.17609... )] 
= 1.1139589... Actual value: 1.1139434... 
Newton’s formula for x = 1 gives 


| ieee 
log 13 = 4 (log 12 + log 14) + 


log 12 = 1.07918... 
log 14 = 1.14613... 
average = 1.11265... 
d = .03347... 
log 13 = 1.11265. . . (6003547 . . . )/26 
= ) ee... 
We find here that Newton’s formula is in fact a superior 


approach to the standard Lagrange interpolation performed on 
twice as many data points! 


An Iterative Relation 


One other approach to finding logarithms involves an iterative 
method known as Borchardt’s Algorithm, accelerated by a tech- 
nique called Richardson extrapolation [8,9]. While this method 
involves extracting a square root and dividing by a multidigit 
number, I find it impossible to exclude because of its extreme 
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accuracy and the minimal memorization required (log 2 and the 
formula itself). The square root is manageable by the techniques 
of Chapter 3; the division may be performed by the cross division 
technique described in Chapter 2 if desired. Perhaps in a pinch 
this provides a good pencil and paper approximation. 

The iterative relation, which we will truncate very early, is 
given for N > 0 and n = 0,1,2,...as 


ap = 4 (1 +N) 
£0 = NI? 
An+1 = Ms (a, mm gn) 
En+1 = (aie . gy)" 


The acceleration is achieved by the relations: 


d(0,n) = a, 
atin) — 2-* @ dik - 1n —- 1) 
~~ 
ee k= 1)2)5,...,n. 

Finally, 

Nem i 
In N = 

d(n,n) 


This seems incomprehensible, but we will only keep terms 


through d(1,1): 
d(0,0) = ap = 4 (1 + N) 


d(0,1) = a, = 4 (ap + go) = 4 (4 (1 +N) +N" 
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d(0,1) — 2 ¢ d(0,0) 


a1). 

(1,1) eee 
N +1 + 4N12 
4 6 


and we end up with the formula, 


OT 19 
nO ON + 1 + 4NIP as 


To use this formula effectively, we need to reduce the range of 
the argument N by extracting powers of 10 and 2 until 1/(2)!” 
< N' < (2)!. This technique of domain reduction is crucial to 
many computer routines using rational or polynomial approxima- 
tions to functions. Since N' = 10°? e 2™ e N, then log N = p + 
m log 2 + (.4343...) ¢ In N'. Of course, the algorithm when 
applied in computers only extracts powers of 2, as these are fast 
bit-shifting operations. 

A plot showing the absolute error in In N' over its reduced 
range is given in Figure 1. The accuracy is astonishing and does 
not require memorization of any constant other than log 2 = 
3010300... (which is easily multiplied by a small integer m) and 
log e = .4342945.... 

Let’s find log 23: 


log 23 = log (10! e 2! e 1.15) 


From the methods given earlier for extracting square roots, we 


find (1.15)!/2 = 1.0723805 .... Then, 


6(.15) 
lo 1. "= 
1.15 + 1 + 4.2895221 


39087 
~ 6.4395221 


= .0606978 ... 
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Figure 1. Absolute error of the logarithmic approximation (Equation 19) 
vs. N. 


log 23 = 1 + .3010300 + .0606978 
= 1.3617278 Actual value: 1.3617278... 


While log 23 is best taken using other methods in this chapter, 
this example serves to illustrate the execution and accuracy of 
the technique. 


Approximate Logarithmic Inverses 


The techniques of the previous sections can produce excellent 
approximations to logarithms of a number, but they provide no 
direct means of finding the inverses, or antilogarithms. Here we 
will consider this inverse as 10*. We should realize that e* can be 
calculated as 10°4343---)*. the latter is much easier to compute. 

A more systematic procedure for finding logarithms can provide 
in tandem a procedure for finding 10%. One method of this sort 
is described by Feynman [10], based on Briggs’ original work. We 
first construct Table 7 for 2th roots of 10. 

Feynman extends the table to the power 1/1024, but | have 
truncated it in an attempt to balance the required memorization 
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Table 7 
2"th Roots of 10 


d = Fractional 
Power Part of 64 109/64 


32 3.16228. 4 
L.(/@Eo. . 
1.339923 
1. 1540S). 28 
1.07461... 
1.08603. .; 


and the desired accuracy. In addition, convenient fractions repre- 
senting the right-hand column involve larger integers as the 
numbers approach 1. The following discussion is valid for any 
length of the table. 

Notice that as the right-hand column progresses toward 1, the 
decimal part of the successive square roots approach that obtained 
by simply halving the decimal parts of the preceding entry. This 
is because, from the binomial series, 


2 
(ee ee 
2 6 


By computing the first 27 roots of 10 to fourteen places, Briggs was 
able by a more involved process described by Goldstine [3] to 
derive the next 27 without loss of precision. 

For small d, we can use the relations 


x #2 


] +i cee tT .. . 
: » 


10* = e (2-303 re) 


to find that 
104 = 1 + 2.3d (20) 
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The factorial “n!” signifies the product of all integers between n 
and 1, inclusive (e.g., 2! = 2 © 1, 3! = 3 e 2 © 1, and so forth). 
The general procedure for calculating log N, then, is to first 
extract powers of 10, say 10°, until N < 10. Then values of 104 
corresponding to dj), d),...are divided out until the absolute 
value of N is less than 10!/4. Then the final term d, can be 
estimated from Equation 20. Finally, since we know that 


N = 10? e@ 1041 + d2 +---+ dnd/64 
then 
log N = p + (d,; + d, +... + d,)/64 


For example, for a = 1211, we extract 10’ to reduce N to 1.211. 
Then, from Table 7, d, = 4 provides the highest 10°/°4 less than 
N, so we divide by 1.15478: 


L2aa 
d, = 4: ——— = 1.04868 
1.15478 
Continuing, 
1.04868 
2=1: = 1.01162 
1.03663 
Now we estimate d;: 
1 —“kOldG2 
d; = ——————- = .00505 
2.30 
Then, 


log 1211 = 3 + (4 + 1)/64 + .00505 
= 3.08318 Actual value: 3.08314... 


Now this seems rather complicated, and it is when compared to 
our earlier methods. The division by 2"th roots of 10 is in practice 
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prohibitive. However, by using the method given in Appendix A, 
we can easily arrive at convenient fractions representing these 
values, easing the calculations. The rational approximations given 
inTable 8 are judgement calls of those giving reasonable accuracy 
and offering small integers with good factorability for simplify- 
ing calculations. 

Now the calculations become a series of simpler multiplications 
and divisions. As it happens, we picked out the two most difficult 
approximations, but the calculation is manageable. We need only 
carry four decimal places through the calculation since the 
rational approximations are not extremely accurate. 


d; = 4: 1.211(84/97) = 101.72/97 = 1.0487 


d, = 1: 1.0487(27/28) = (1.0487/28) © 9 @ 3 = 1.0112 


d, = .0112/2.30 = .0049 


log 1211 = 3 + (4 + 1)/64 + .0049 = 3.0830 


which suffers some from the approximations. 

At any rate, the procedure is complicated for finding logarithms 
and requires some memorization. I present it solely because the 
inverse of this method provides a means of finding 10*, which is 
generally a more difficult task than finding log x. 


Table 8 
Rational Approximations to 2"th Roots of 10 


Rational 
Approximation 


3.16228 253/80 
1.77828 16/9 

1.35322 4/3 

1.15478 97/84 = 97/(3 « 4 © 7) 
1.07461 43740 

1.03663 28/21 = 4 « 7)/GP 
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To illustrate this, let us reverse the problem and calculate 10* 
for x = 3.08314. First, we extract the integer portion: 


1093-08314 ‘a 10° e 10:08314 


Then we extract inverse powers of 2 (1/2, 1/4, 1/8, etc.) until we 
reach a value less than 1/64. Unless these fractions are memorized 
or calculated on the run, it may be best if possible to multiply the 
argument by 64: 


(.08314)(64) = 5.3210 
Now 5 = 4 + 1 and 1072! ~ 1 + 2.30(.32)/64, giving 


103-98314 = 103 e (10'4+1)/64) @ 10:3210/64 


u 


10° © (97/84) © (28/27) © (1.0115) 
= 10° e (97/81) @ (1.0115) 
2 1211 to four digits 


Another version of this general method was given by Bemer in 
1958 for use in early decimal computers [11]. The idea here for 
finding logarithms is to first extract powers of 10 until N is less 
than one. Then, N is multiplied by a series of simple numbers 
until N is just greater (or less) than 1, at which point the formula 
for In (1 + x) given in Equation 14 is invoked. Bemer uses the 
multipliers shown in Table 9, based on the first digit following the 
decimal point in the current value of N (other multipliers may be 
found in the work of Camp [12]). 

Here the multipliers are simply obtained as the highest two- 
digit number that, when multiplied by the maximum value in the 
range, will result in a new value of N < 1.1. For example, the range 
from .1 to .2 would have a maximum value of .2, so the multiplier 
is 1.1/.2 = 5.5. A value of N that will lie between 1.0 and 1.1 
will be reached in at most three multiplications. The series for 
In (1 + x) is then taken to the desired accuracy. 

For our purposes, it is probably easier to directly use the 
multipliers 1.1/.2, etc., as shown in Table 10. 
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Table 9 
Multipliers in Bemer’s Method for Logarithms 


N to One Digit Multiplier 


0 MNDAWNBWN Ye 


Table 10 
Multipliers in Our Method for Logarithms 


N to One Digit Multiplier Log (Multiplier) 


log 11 — log 2 
log 11 — log 3 
log 11 — log 4 
log 11 — log 5 
log 11 — log 6 
log 11 — log 7 
log 11 — log 8 
log 11 -— log9 = 
log 11 — log 10 = 


l 
vE 
5, 
7 
fe 
fe 
ar 
8 
iv 


Again, we find log 1211: 
log 1211 = 4 + log .1211 
.1211(11/2) = .66605 

.66605(11/7) = 1.0466 
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log 1211 = 4 — .7404 — .1963 + .43(.0466) 
= 3.0833 


Of course, greater accuracy can be obtained by using more digits 
in the logarithms and more accuracy in the series expansion. 
Actually, we could have noticed that the intermediate value 
.66605 of N can be brought very close indeed to 1 if multiplied 
by 1.5, a logarithm easily obtained from those of low primes. It 
pays to be flexible. 

Finding 10* with this method parallels the method using even 
roots of 10. Powers of 10 are extracted from the number N until 
N < 1. Then the logarithms of the multipliers are subtracted until 
N is very near 1. If we do not care to explicitly memorize these 
logarithms as unique numbers, we can in each step subtract log 11 
= 1.04139 and add the logarithm of a single-digit number (which 
we should be able to quickly obtain), leaving the result near 1. 

The procedure continues as before: 


1032-08314 = 103 e 10-98314 
= 10° ¢ (11/9) « 10°! 
From Equation 20, 
10-4! = 1 + 2.3(-.0041) = .9906 
giving, to four digits, 
10308314 ~ 121] 


Now let’s try another example, 10!°’7!° = 47: 


10167210 =l0l e 1067210 


10! © (11/3) « 10-1078 


10! © (11/3) © (11/9) « 10° (21) 
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Again, 

10-0206 ~ 1 + 2.3(.0206) = 1.0474 
and we arrive at 

10167210 - 46.94 

Once again, greater accuracy is possible with more work. Given 
that log (11/10) = .0414, which can actually be added or sub- 
tracted, we ended up within 1 part in 10,000 of the worst (i.e., 
maximum) power of 10 we could have (+.0207). Therefore, we 
know we can only improve our accuracy if we choose, say, 11/4 


as the first multiplier instead of the prescribed 11/3: 


101-67210 = 10! e 109-6721 
= 10! ¢ (11/4) © 10238 
= 10! © (11/4) © (11/7) © 10°36 
= 10! © (11/4) © (11/7) © 10-%9 
= 10! © (11/4) © (11/7) © (.9887) 
= 47.00 


We can work out situations where the exponent of 10 ends up 
relatively large, say of a magnitude greater than .005. These cases 
amount to 10'@*5), where a = .02, .01, —.01, or —.02 and b < .005. 


Consider the case a = .02: 
10-2 = 1.04713 
1 + 2.3(.02) = 1.04600 


difference = .00113 
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Then, 
1 Q{-02+b) . 10-°2 e 10° 
= {1 + 2.3(.02) + .00113] © (1 + 2.3b) 
= 1+ 2.3(.02 + b) + .1b + .000113 


Also, 
10-(-02+b) = 1 — 2.3(.02 + b) + .1b + .00099 


Therefore, we will significantly improve our approximation of 
10*@*) if we always add to the prescribed quantity 1 + 2.3(a + b) 
the term (cb + d), where c and d are given in Table 11. 

This allows us to achieve greater accuracy from our earlier 
intermediate equation (Equation 21) without much additional 
work. Now while we could memorize 10°’, ... , 10~° and directly 
extract these multipliers as needed, I think this is a more difficult 
task of memorization. More importantly, this would leave us with 
two multidigit numbers to multiply in the final result (10° and 
10°) rather than one. 

We arrived earlier at Equation 21: 


101-67210 = 10! e (11/3) © (11/9) @ 10020 
We now have: 
10-9206 = 1 + 2.3(.0206) + .1(.0006) + .00113 
= 1.0486 


Table 11 
Correction Values for 10:!¢+>) 


OOO29 
002 4 
.00099 
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10!-67210 ~ 10! @ (11/3) © (11/9) © (1.0486) 
= 46.99 


We were actually limited by the four-place precision we retained 
throughout the calculation. The example given in the next 
section maintains five digits throughout. 


An Example Problem 


To demonstrate the use of logarithms in extracting higher-order 
roots, let’s return to the mental calculation mentioned at the 
beginning of this chapter, the seventh root of 13781. I picked the 
number 13781 from a table of prime numbers and made every 
effort to avoid one offering special qualities for simplifying the 
logarithmic conversions. This is actually somewhat of an effort as 
the overwhelming majority of numbers have some sort of unique 
quality, such as being located near a number that is easily factored. 

At any rate, we must first find log 13781. Personal bias leads me 
to the following approach: 


log 13781 = log 14000 + log (1 — 219/14000) 
log 14000 = 3 + log 2 + log 7 
= 3 + 30103 + .84510 
= 4.14613 
Because 219/14000 is not really very small, we need to retain two 


terms in Equation 14, although the second can be calculated to 
very low precision. 


1 Zi9 LY 2s 
log (1 — 219/14000) = — - = = [| 
2.3L 14000 2 \14000 


219/14000 = .01564 
¥ (01564)? = ¥% (15)(16) « 10% = 12 e 105 


es 
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log (1 — 219/14000) = (-.01576)/2.3 = —.00685 
Therefore, 

log 13781 = 4.13928 

(1/7) log 13781 = .59133 

Notice that the order of the root is almost inconsequential here. 
We now need the antilogarithm of this result. Using Bemer’s 
method and retaining logarithms of the multipliers to five places, 
we find: 

log 11 — log 3 = .56427 

1059133 = (11/3) © 1002706 

log 11 — log 10 = .04139 

10:59133 = (11/3) © (11/10) ¢ 10--91433 


Now, from Table 11, which we memorized if we desire extreme 
accuracy, 


1091433 ~ 1 + 2.3(-.01433) + (.1/2)(.00433) + .00024 
= .96750 
which yields 
1079133 = (11/3) © (11/10) © (.96750) 
= 3.90225 
compared to the actual seventh root of 13781, 3.90235.... 
There are many other ways, of course, to determine log 13781; 
upon reflection, I think a number of them will become apparent. 


Again, most numbers can be calculated without the second term 
in Equation 14 and possibly without the corrections given in 


Table 11. 
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Since the number of significant digits of a logarithm is the same 
as that of the number whose logarithm is taken, we need only 
consider five digits of the number if a five-digit result is desired. 
The value log 1378131, for example, may be represented as (2 + 
log 13781) to this accuracy. 


Bibliography 


1. E. G. Kogbetliantz, Fundamentals of Mathematics from an 
Advanced Viewpoint Vol. II: Algebra and Analysis: Determi- 
nants—Equations—Logarithms—Limits, Gordon and Breach, New 
York, 1968, pp. 470-490. 

2. Steven B. Smith, The Great Mental Calculators, Columbia 
University Press, New York, 1983, pp. 150-155. 

3. Herman H. Goldstine, A History of Numerical Analysis From 
the 16th Through the 19th Century, Springer-Verlag, New York, 
1977, pp. 1-62, 70-71. 

4. N. David Mermin, “Logarithms!,” American Journal of Physics, 
46 (1978) pps IDI- 15, 

5. William R. Ransom, “Elementary Calculation of Logarithms,” 
The Mathematics Teacher, 47 (1954) pp. 115-116. 

6. Albert A. Bennett, “Note on the Computation of Logarithms,” 
American Mathematical Monthly, 28 (1921) pp. 130-131. 

7. James C. Kirby, “An Efficient Logarithm Algorithm for Calcu- 
lators,” College Mathematics Journal, 19 (1988) pp. 257-260. 

8. B. C. Carlson, “An Algorithm for Computing Logarithms 
and Arctangents,” Mathematics of Computation, 26 (1972) 
pp. 543-549. 

9. George Miel, “Of Calculations Past and Present: The Archi- 
medean Algorithm,” American Mathematical Monthly, 90 
(1983) pp. 17-35. 

10. Richard P. Feynman, The Feynman Lectures in Physics Vol. I, 
Addison-Wesley, Reading, 1963, section 22-4. 

11. R. W. Bemer, “A Subroutine Method for Calculating Loga- 
rithms,” Communications of the Association for Computing 
Machinery, 1 (1958) pp. 5-7. 

12. C. C. Camp, “Logarithms of Large Numbers,” American Math- 
ematical Monthly, 35 (1928) pp. 547-551. 


‘Trigonometric 
Functions and 
Their Inverses 


Although this may seem a paradox, all exact science 
is dominated by the idea of approximation. 


Bertrand Russell [1] 


—- functions offer unique challenges to us in our 
attempts to provide convenient techniques for mental calculation. 
For one thing, their derivatives yield other trigonometric func- 
tions, so efforts to use approximations based on derivatives (such 
as the Newton-Raphson method) are generally fruitless. The 
orthogonality of the sine and cosine functions, which allows us to 
represent any function as a series of sine and cosine terms, 
correspondingly makes it difficult to represent them with conven- 
ient nontrigonometric functions. In addition, arithmetic opera- 
tions such as sin (a + b) or (sin a + sin b) do not lend themselves 
well to simplification. Further, the values of at least the sine and 
cosine functions lie in the range from —1 to 1, sharply reducing 
the validity of deleting higher-order divisions or truncating 
power series. The tangent function spans the other extreme, from 
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—oo to +co. Finally, the radian unit of angle measure (in which the 
functions are most naturally and simply expressed in relations) 
spans for the first O° to 90° quadrant an almost equally maddening 
range of 0 to 1.57. 

Despite these difficulties, we will discuss here some methods of 
use in reckoning the basic trigonometric functions sin x, Cos x, 
and tan x, as well as their inverses arcsin x, arccos x, and arctan x. 
Only the first quadrant is considered, as trigonometric values in 
other quadrants are easily deduced as sign changes from those in 
the first quadrant. 


Sine and Cosine Functions 


We may approach the sine function by utilizing the familiar 
power series: 


e 0 


sm x = x-— +] —-. 


3! #8 


where again the factorial function n! represents the product of all 
integers between n and 1, inclusive. 

We can truncate this series, using “relaxed” coefficients to yield 
less overall error in this formula [2]: 


sin x ~ .99989x — .16596x? + .00760x? 


Now suppose we know the sine of an angle a in degrees and we 
wish to find the sine of an angle d a small number of degrees b 
from a. We can then rewrite the last relation for sin (a + b), where 
the radian equivalent x to (a + b) is given by n(a + b)/180. The 
Greek letter pi is the famous constant of proportionality between 
the circumference and the diameter of a circle, and equals 
3.14159265 .... 1 assume throughout this chapter that we work 
in degrees instead of radians, as this occurs an overwhelming 
portion of the time in practice. If a calculation in radians does 
arise, we can convert the argument to degrees by multiplying it 
by 180/z, a quantity conveniently approximated by the fraction 
401/7, acquired by the means detailed in Appendix A. 
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To continue, we arrive at 
sin (a + b) = .017451(a + b) — 8.8234 © 10°-"(a + bp 
+ 1.23 e 10°'!(a + b)? 


We can multiply this out, extract the approximation for sin a 
from the right side and drop higher-order terms in b as negligible. 
Multiplying through by 1000 to ease calculations then yields 
the approximation, 


b 
1000 sin d = 1000 sin a + — (174 _ ~) (22) 
10 40 


where d = a + b. 

This formula is arranged to simplify the mental computation 
and is actually far easier to calculate than it first appears. It has 
been very useful to me, and it represents my best attempt at a 
reasonably convenient approximation. 

How good is it? The jagged curve plotted in Figure 2 represents 
the absolute error curve, or the quantity (approximation — func- 
tion) for this relation over the first quadrant. This curve assumes 
that sin a is known for all angles that are a multiple of 10°, 
yielding —5° < b < 5°. 

Looking at the jagged curve, we can see that for angles d 
through 54°, represented by a = 50° and b = 4°, the approximation 
is accurate to three (and usually four) decimal places. Theretore, 
we can use this formula in the range 0° < d < 54° if we memorize 
the following values: 


sin 0° =0 
sin 10° = .1736 
sin 20° = .3420 
sin 30° = 5 
sin 40° = .6428 
sin 50° = .7660 


148 = Dead Reckoning Calculating Without Instruments 


0 10 20 30 40 50 60 70 80 90 


Figure 2. Error curves (approximation — function) for Equation 22, given 
by the jagged plot, and Equation 23, the smooth plot, vs. angle in degrees. 


A useful aid in jogging the memory is to recall from Equation 22 
that the initial slope of the sine function is 0.174. The value of 
sin 10° is just slightly lower than this because the curve begins to 
flatten. The value for sin 20° shows some additional flattening. 

For those who do not wish to memorize any values of sin a at 
all, we can now backtrack and find a very reasonable approxima- 
tion of the sine function in the same interval 0° to 54°, although 
the average error will be significantly greater than that just given. 

In the derivation of Equation 22, if we subtract a quantity that 
is relatively small over the range given, namely 


b (= + ") 
10 \ 720 
we will find that we can derive the approximation, 


d(d + =| 
120 


d 
1000 sin d = 7 [1744 = (23) 
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where we have now kept four places in 174.4 because it is more 
significant here than in Equation 22, where the known values 
sin a “pulled” the curve closer to the correct one. 

The error curve for this approximation is also given in Figure 2 
as the smooth plot. The error is quite low over the range 0°—54°, 
but nonetheless it is greater on average than the other. 

It is important to recognize that the plots in Figure 2 are 
designed to accentuate the errors. If the true sine function from 
O° to 54° were plotted on a full page here, the difference between 
it and either of the sine approximations (Equations 22 and 23) 
would be less than the width of the printed line itself! 

Now returning to the quandary that began this entire enter- 
prise, as mentioned in Chapter 1, we can easily find sin 28°. From 


Equation 22, 


30 @ 28 
1000 sin 28° = 1000 sin 30° — .2 (: 74 - ——7) 
= 580 — 30:6 
or, 
sin 28° = .4694 Actual value: .46947... 


From Equation 23, 
a3 ° 79 
1000 sin 28° = 2.8 (174.4 = eS 
120 


sin 28° = .46937 


Let us defer the task of finding sin d for d > 54° and instead 
address the cosine function. We begin with the identity, 


sin (a + b) = sin a cos b + cos a sin b 
or, 


sin (a + b) - sin a cos b 
Cos 3 ee 
sin b 
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For b = 1°, then sin b = .01745 and cos b = .99985. Substitut- 
ing our approximation (Equation 22) for sin (a + b), and setting 


b = 1° yields 
cos a = 1 — .00014327a(a + 1) 


Changing a to the variable d and realizing that 1/7 = .1429..., 
we find 


d(d + 1) 
1000 cos d = 1000 — > (24) 


The smooth curve in Figure 3 displays the absolute error (again, 
approximation — function) for this relation as a function of angle. 
This approximation shows reasonable accuracy (usually three 
digits, but rarely more than one off in the third digit) through 40°. 
A passing familiarity with the hump at around 30° can result in 
three-digit accuracy throughout this range. This is the formula | 


0 10 20 30 40 50 60 70 80 90 


Figure 3. Error curves (approximation — function) for Equation 24, given 
by the smooth plot, and Equation 25, the jagged plot, vs. angle in degrees. 
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have used on occasion, and it requires no memorization of cosine 
values. Using this formula without recalling that we are in the 
hump area, we find that cos 28° is approximately .8840, compared 
to the actual value of .8829.... 

We can take this formula, however, and form a more uniformly 
accurate approximation if we memorize the cosines of angles that 
are again multiples of 10°. For d = a + b, Equation 24 becomes 


(a + b) @ (a + 1 + b) 
7 


1000 cos d = 1000 — 


Substituting Equation 24 again for cos a, 


b(a + d+ 1) 
1000 cos d = 1000 cos a - > a (25) 
The jagged trace in Figure 3 comprises the error curve for this 
function. It shows an improvement on average over that shown 
for Equation 24 in the interval 0° < d < 35°, assuming that for 
d = 35° we take a = 30° and b = 5°. To use this relation, we 
therefore need to know the following values: 


cos O° =1 
cos 10° = .9848 
cos 20° = .9397 


cos 30° = V¥3/2 = .8660 
With this relation, we find cos 28° as: 


(2) @ 59 
1000 cos 28° = 1000 cos 30° —- ——————_ 


cos 28° = .8829 


which is accurate to four digits. 
Nature has been unusually kind to us in both of these relations 
because the range of angles for which the sine approximations 
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(Equations 22 and 23) are invalid is almost exactly (or is exactly) 
the range of angles for which the cosine approximations (Equa- 
tions 24 and 25) are valid, and vice-versa. Therefore, we can let 
sin d for d > 54° become cos (90° — d), whose argument is less 
than 36°, and cos d for d > 35° become sin (90° — d). In the end, 
we cover the entire 0° to 90° range of interest for both the sine 
and cosine functions with two formulas. Incidentally, the alter- 
native if this were not true would be to add a third relation, an 
identity, to double the useful range of the cosine approximation: 


cos 2d = 2 cos? d—-1 


The corresponding double-angle identity for the sine function 
mixes sine and cosine terms and is not useful. 

We may note in closing this sine-cosine discussion that Equa- 
tions 23 and 24 yield the approximation, 


172 sin d 
= ———_——__ - 2 
d 


COS 


For a right triangle with sides of length a < b < c, we can formulate 
the above relation to express the smallest angle A of the triangle 
in terms of the three sides [3,4]: 


a 


A = 172 (26) 


+ 2c 


where A is in degrees. The third angle is then apparent. 

This approximation was given by Ozanam in 1699 and is 
actually very good, giving almost four-digit accuracy for angles up 
to 45° (where the angle is no longer the smallest). For the familiar 
3,4,5 right triangle with a smallest angle of 36.870°, this 
formula gives 


A = 172 e (3/14) = 36.857° 


The Tangent Function 


The tangent function, of course, may by definition be calculated 
by dividing the sine function by the cosine function for the given 
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angle. Since these formulas are not valid in the same range, 
though, the calculations reduce in practice to: 


sin d 
tan d= tar 0° < d < 36° 
cos d 
sin d 

= for 36° <d < 54° 
sin (90° — d) 
cos (90° — d) 

= —___—— for 54° <d < 90° 
sim (90 - d) 


These need not be remembered, as the proper terms become 
apparent when performing the individual sine and cosine 
calculations. 

Using the sine relation (Equation 22) and cosine relation (Equa- 
tion 24), we find tan 28° to be .5310... instead of the actual 
value of .53171.... Using the cosine relation (Equation 25), we 
arrive at a value of .53168.... 

There exists a power series for tan x that can be truncated with 
relaxed coefficients [5]: 


tan x = x + .31755x? + .20330x° 


The magnitude of the error in this approximation is less than 
0.001 for the range 0° < x < 7/4, or 45°. However, the relatively 
large coefficient of the last term indicates difficulty in further 
simplification. If we proceed in the same manner as for sin x, 
though, recognizing that the form of the relation is the same, we 
derive for d = a + b, 


b d 
1000 tan d = 1000 tan a + — (174 + ~) (27) 
10 20 


where we memorize or recall the following values: 
tan 0° = 0 


tan 10° = .1763 
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tan 20° = .3640 

tan 30° = V3/3 = 5774 
tan 40° = .8391 

tan 45° = 1 


Note the change in sign of the last term in Equation 27 compared 
to the earlier sine approximation (Equation 23) of the same form. 
The error curve for the approximation (Equation 27) is shown 
in Figure 4. The increase in difficulty in approximating tan d 
compared to sin d or cos d is apparent. 
Equation 27 can be modified empirically to greatly improve the 
accuracy, but with additional memorization required: 


nb? 


fee 
40 


b d 
1000 tan d = 1000 tan a + — (174 mn + _ (28) 
10 20 


0 5 10 i> 20 Zp 30 35 40 45 


Figure 4. Error curve (approximation — function) for Equation 27 vs. angle 
in degrees. 
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The values of m and n appear below, where they have been 
contrived to aid memorization: 


gl 0 
wm be A 
a SC! 
. 15 3 
40 45 10 


Figure 5 displays the error curve for this approximation, 
which is for all practical purposes accurate to at least three 
decimal places (or + 0.0005) for 0° < d < 45°, and requires no 


significant division. 


35 40 45 


0 5 10 15 20 25 30 


Figure 5. Error curve (approximation — function) for Equation 28 vs. angle 


in degrees. 
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Using this formula, 


a 304 
+ 


1000 tan 28° = 1000 tan 30° + (-.2) © (189 4 40 


tan 28° = .5504a. .. 


compared to the actual value of .53171.... 

The range 45° < d < 90° is very difficult to approximate 
explicitly because the tangent function tends to infinity. However, 
for this range we can use the trigonometric identity, 


1 


tan d = ———————_ 
tan (90° — d) 

and suffer the division, although with this disadvantage it is 
probably easier to take the sine-cosine ratio. 


The Arcsine and Arccosine Functions 


The trigonometric inverses are exceedingly difficult to approx- 
imate well mentally, but are not often required in practical work. 
Looking at Equations 22 through 28, the only one amenable to 
inversion is the simpler cosine function (Equation 24). The result- 
ing expression for arccos x after some simplification is given by 


arccos x = [7(1000 — 1000x)]!/2 — .5 (29) 


Three identities linking the arcsine and arccosine functions are 
of use to us here: 


arccos x = arcsin (1 — x’)!/ (30) 
arcsin x = arccos (1 — x’)!/2 (31) 
arcsin x = 90° — areg@es x (32) 


We can then substitute the last two of these into Equation 29 to 
arrive at another relation for arccos x: 
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arccos x = 90° — [7(1000 — 1000(1 — x*)!/2)]!/2 _ 5 (33) 


Plots of the absolute error curves for Equations 29 and 33 are 
given in Figure 6 for O < x < 1, or 90° > d => 0°, where d = arccos 
x in degrees. As is evident, Equation 29 provides best results for 
x above .707 (d < 45°) and Equation 33 for the complement. The 
two relations coincide at d = 45°, where x = (1 — x?)!/. 

This pair of formulas yields angles accurate to + 0.5° over the 
first quadrant (other quadrants are easily deduced). The square 
roots are readily found to four or more digits by the methods 
described in Chapter 3. In fact, the squaring operation here should 
be more intimidating at this point than the square root extraction. 
The arcsine can also be easily determined through these relations 
and the Equations 30-32. 

Can we use our memorized values of the sine function to aid 
in calculating the arcsine (and arccosine) functions? Given that 
the reader is interested in this abstruse section, we may assume 
that sine values have been memorized from the last section for the 


0 0.25 0.5 0.75 1 


Figure 6. Error curves (approximation — function) for Equation 29, the 
rightmost plot, and Equation 33, the leftmost plot, vs. the cosine of angles 
to 90°. 
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first several multiples of 10°. Now, a series relation for arcsin x, 
where —1 < x < 1, exists in a similar form to those of the last 
section [5]: 

x? iwsen 1-3 5 om 


2 re 
i Ee Swe. ae 


arcsin X = x + 


Truncating this series after the second term, substituting x = a + b, 
dropping a b* term and separating out the terms comprising 
arcsin a, we find 


ax 


401 
arcsin (a + b) = arcsin a + = b (1 + =) (34) 


The familiar coefficient 401/7 converts the result into degrees. 

The general procedure, then, is to recall the sin d value nearest 
x and call it a. Then b = x — a. As in general a is not a round 
number, and because we multiply by 401/7, Equation 34 is 
somewhat more difficult than the corresponding formulas from the 
last section. 

If x = .46947, then a = .5 = arcsin 30° and b = —.03053. 
Then, 


01 5)(.469 
acsin x = 30° + © (0305) « [1 « PACH 


2 
== 30° — 1.95° 
= 28.05° 


compared to the actual value of 28°. 

The error curve for Equation 34 is shown in Figure 7. As 
expected, the approximation becomes poorer for larger x, even 
though the intervals between values of a become smaller. The first 
45° range (0 < x < .707) is quite acceptable, however, and 
Equations 30-32 then provide arcsine and arccosine values 
throughout the first quadrant. 
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Figure 7. Error curve (approximation — function) for Equation 34 vs. the 
sine of angles to 90°. 


For additional accuracy, we find from Hastings [2] that 
arcsin x = 1/2 — (1 — x)! @ (1.57073 — .21211x 
+ .07426x? — .01873x?) 


with an error no greater than 0.00005 for 0 < x < 1. 

The second term is obviously that for arccos x. Since we desire 
x to be as small as possible, let us consider as before the range 
O <x <.707. We can truncate this expression, realizing that we 
need to multiply by 180/ to convert arccos x into degrees d. If 
we drop terms involving powers of x greater than 1, and empiri- 
cally add a correction term on the end of our resulting relation, 
we find 


arccos x = (1 — x)!/2 e (90 — 12x) + 1.5(x —.1) (35) 


A plot of the error function for this approximation is shown 
in Figure 8. We deduced earlier from Equations 31-32 that 
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Figure 8. Error curves (approximation — function) for Equation 35, opti- 
mized for x < .707, and after modifications, optimized for x > .707, vs. the 
cosine of angles to 90°. 
arccos x = 90° — arccos(1 — x’)!. The corresponding error 
function for this approximation is shown as well in Figure 8, where 
Equation 35 is again used for approximating the arccosine func- 
tion. This latter approximation represents the arccosine function 
with greater accuracy for x 2 ./07. 

With Equations 30-32, we have covered the first quadrant for 
the arcsine and arccosine functions to reasonable accuracy, 
given the infrequency of the situation. Taking again the example 


x = .4694716, 
arccos x = 62.0044° 


compared to the actual value of 90° — 28° = 62°. 

It is also worth mentioning more exact approximations (in some 
regions) for arcsin x and arccos x formulated from our earlier 
approximation (Equation 26) for the smallest angle A in a right 
triangle of sides a < b < c. If we define a = x and c = 1, then 
A = arcsin (a/c) = arcsin x. Writing b = (c? — a2)!/? = (1 — x*)!? 
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and replacing the number 172 by 3(401/7), we arrive at 
the expression, 


5x 401 
arcsin X = reo e (—) (36) 


Likewise, setting b = x, c = 1, anda = (1 — x’)!/", we find 


(3#) 


arccos X = 


ae ) 
ae 2 | 

An alternative means of deriving these relations uses the 
accelerated Borchardt’s Algorithm described in Chapter 4 [6]. To 
the iterative relations given there, we apply the initial values 
ay = (1 — x’)! and g, = 1 and produce approximations for the 
arcsine function, the first being Equation 36. Setting ay = x and 
go = 1 produces approximations for the arccosine function, the 
first given by Equation 37. 

We find here that the result for the arcsine function is identical 
to that obtained by transforming x in the arccosine formula by 
(1 — x’)!”. The error curves for Equations 36 and 37 are plotted 
in Figure 9. Choosing the more accurate equation and, if neces- 
sary, using the earlier identities gives the best accuracy yet for 
these functions, albeit at the expense of more difficult computa- 
tions. Perhaps these relations are best left to cases where x is a 
value of one or two digits. 


The Arctangent Function 


One approach to approximating the arctangent function utilizes 
the following relation [5], with the multiplier again added for 
conversion to degrees: 


401 
arctan x = —— ® oe (38) 


ij 1 + .28x2 


for-1 <x <1. A plot of the corresponding error function is shown 
in Figure 10, showing an absolute error less than 0.28° over the 
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Figure 9. Error curves (approximation — function) for Equation 36, the 
leftmost plot, and Equation 37, the rightmost plot, vs. x. 


0 Own25 On.5 0.75 L 


Figure 10. Error curve (approximation — function) for Equation 38 vs. the 
tangent of angles to 45°. 
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range —1 < x < 1, or O° <d < 45°, where d = arctan x in degrees. 
A potential difficulty exists for the range 45° < d < 90°, however, 
as the arctangent function careens to infinity. The solution is to 
use the identity, 


arctan x = 90° — arctan (1/x) (39) 


to reduce x in this region to a value less than 1. To avoid excess 
computations, we can replace x in Equation 38 with 1/x and 
simplify, giving 


op 3 _* (40) 
arctan xX = a 
( war + 28 


iat & > I. 

Another method makes use of memorized values of tan d for 
d = 0°, 10°, 20°, 30°, 40°, and 45°, as given in the discussion on 
the tangent function. If these values are memorized for calculating 
tan d, they may be used here as well to great advantage. 

The derivation follows from an identity for the tangent function 


[/). For & = a+ b, 


tan d = ® 
tane=a 
tan d — tan e 7a 
tan (d — e) = ———————_ = 
1 +tandtane 1 + ax 
or, 


b 
d = arctan a + arctan [——_ 
1 + a(a + b) 


For small arguments of the last term, we arrive at 


401 b 
arctan x = arctan a + —— ® 
7 1 + ax 


(41) 


forx =a+b. 
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A plot of the error function for 0 < x < 1 is given in Figure 11. 
Notice the order of magnitude change in the scale; this is an 
extremely accurate approximation with little increase in complex- 
ity over the last, given the memorization required. Again, Equa- 
tion 39 is used for x > 1. 

For those interested, Equation 26 can be used to derive an 
extremely accurate approximation to the arctangent function as 
well as the earlier arcsine and arccosine functions. We derive: 


3x 401 
Eee ee cor = 
1+20+x)”" 7 


arctan x = 
The accelerated Borchardt’s Algorithm can also be used to derive 
this; the initial values used here for the arctangent function are 
given as ap = 1, g = (1 + x’)!”. 

The extreme accuracy of this relation for O < x < 1 seems 
unnecessary considering the accuracy available from the pre- 
vious relation. 


8) 0.25 0.5 0. 75 i 


Figure 11. Error curve (approximation — function) for Equation 41 vs. the 
tangent of angles to 45°. 
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Other Approximations 


Many other approximations exist for trigonometric and inverse 
trigonometric functions. However, the vast majority do not lend 
themselves even as well as the ones given here to mental (or back 
of the envelope) solution. For those interested, references for 
other approximations include Fike [8], Kogbetliantz [9-12], Fox 
[13], Frame [14], and Spielberg [15]. 

Interestingly enough, mental calculation of hyperbolic and 
inverse hyperbolic trigonometric functions (an extremely impres- 
sive talent) is much easier, since they are simply defined in terms 
of real exponential and logarithmic functions more readily eval- 
uated by the methods of the last chapter. For example, 


2 arctanh x = log (1 + x) — log (1 — x) 


Other relations are easily found from definitions of these functions. 
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Concluding Remarks 


The object of pure Physic is the unfolding of the 
laws of the intelligible world; the object of pure 
Mathematic that of unfolding the laws of 
human intelligence. 


J. J. Sylvester [1] 


To think the thinkable—that is the 
mathematician’s aim. 


C. J. Keyser (1904) [1] 


T then, represents a collection of algorithms specifically 
developed for mentally calculating or approximating arithmetic 
results and elementary functions. I hope that a sense of adventure 
and marvel has been conveyed, and a certain eagerness inspired 
to try out the techniques as the opportunity arises or is created. 

I also strongly recommend looking through the referenced 
articles and books for more detailed information on any topic that 
piques your interest. I appreciate that this may be an untamiliar 
step for some people; it shouldn’t be, and the majority of the 
references can be found in the library of even a small college or 
obtained through any public library. The many articles and hooks 
listed are clear ones and are certainly within the grasp of anyone 
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who appreciates the topics in this book. Some of them, including 
those by Mermin [2], Goldstine [3], Taylor [4], Menninger [5], 
Aitken [6], Loweke [7], Uspensky [8], Eynden [9], Bailey [10], and 
Kovach [11] (in no particular order) are simply a pleasure to read. 
As you may have gathered, Smith’s book on mental calculators 
[12], which has a somewhat different emphasis, complements this 
book very nicely. 

One important fact implicit in the presentation of the methods 
in this book is that they are in no way the final word on the 
subject. I encourage the reader to find better and alternate 
methods. Experimenting with perceived properties of numbers or 
functions can be done at any odd time on any scrap of paper, and 
this is an area where playing with particular numbers often leads 
to general algebraic relationships or approximations. It invariably 
leads to some insight into the intricacies of the number world. 
Hoist a sail and catch the wind. 
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Appendix 


Finding Rational 
Approximations to 
Precomputed 
Constants 


The employment of mathematical symbols is perfectly 
natural when the relations between magnitudes are 
under discussion; and even if they are not rigorously 
necessary, it would hardly be reasonable to reject 
them, because they are not equally familiar to all 
readers and because they have sometimes been 
wrongly used, if they are able to facilitate the exposition 
of problems, to render it more concise, to open the 
way to more extended developments, and to avoid the 
digressions of vague argumentation. 


A. Cournot (1897) [1] 


Cen we encounter multidigit numbers, such as scaling 
factors, which are difficult to use directly as a divisor or multiplier. 
In addition, there are certain types of calculations, such as those 
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for finding antilogarithms, that benefit from range reduction 
techniques, i.e., the original number is reduced in some manner 
into a range where we can obtain a good approximation. We then 
adjust this approximation to find that for the original number. 
This whole process generally involves multiplications or divisions 
by constants (such as 107", where n is an integer) that are 
precomputed and memorized. This appendix describes the general 
procedure by which we can generate whole-number fractions that 
approximate multidigit numbers to a required accuracy and in 
practice involve multiplication and division by reasonably small 
whole numbers. 

To illustrate, to convert from radian units to degrees, we 
multiply by 180/m = 57.296.... This is difficult to do, and in 
Chapter 5 I use instead the multiplier 401/7 = 57.286.... This 
is a fraction containing numbers that are easy to multiply and 
divide by, yet is very close to the actual conversion factor. More 
accurate fractions are 974/17 and 4068/71; worse ones include 
172/3 and 229/4. These approximations are straightforward, but 
tedious, to generate, so I use a computer program to create them. 
I then look over them to find one that is a reasonable approxima- 
tion, but tnat has a numerator and denominator that I can 
mentally multiply and divide by relatively easily. Once a fraction 
(such as 401/7) is selected, I always use it for this constant in any 
calculation that requires it. 

We approach this subject in terms of interesting creatures 
called continued fractions. A continued fraction is a fraction of 
the form: 


ay 


apt 


1 


ager... 


where the integers a, are termed partial quotients. A more modern 
notation for continued fractions is given by the equivalent form, 


1 1 1 


a a ta 


Aye 
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A comprehensive treatment of continued fractions may be found 
in the work of Olds [2]; here we are interested in their application 
to rational approximations [3,4]. 

As examples, it is commonly given that: 


of 9 -_ See 
1+2+1+1+4 + 
1 l l 1 
=3+— — ... 
PS Mah s.2092 + 


For integer values of a, the continued fraction will converge to 
a limiting value. Terminating the fraction at given a,’s and 
simplifying the fraction produces rational approximations to the 
value called the principal convergents A,. 

We need, therefore, a procedure for converting a given number 
into a sequence of numbers a,, that comprise a continued fraction 
for the number. We can extract the principal convergents A,, until 
we arrive at fractions involving numbers too large for us to work 
with. Convenient fractions can then be tested for accuracy to the 
original number. As a rule, we expect good approximations for a 
given complexity of the fraction to occur when the next a, 
following termination is relatively large, as its reciprocal would 
provide a correspondingly small correction if included. For exam- 
ple, in the continued fraction for m given earlier, terminating the 
fraction at a, (before as = 292) gives 355/113, a value accurate to 
7 digits. Finally, we can find intermediate convergents; these occur 
for any principal convergent A,, that results from truncating the 
continued fraction at a partial quotient ap greater than 1. 

It can be shown that the values a, are simply the multipliers 
in each step of Euclid’s Algorithm for finding the greatest common 
divisor of two integers (see Chapter 2). To prove this, consider two 
numbers a and b with a > b. Euclid’s Algorithm gives: 


a=ab+t+r, 
b = apr; + 
Ty > agig * ty 


ete. 
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Dividing each line by the variable associated with the a, term, 
we find 


a Tj 
b b 
b 12 
T} rT} 
Tj T3 
19) r2 
ete; 


Then, since the last term in each line is the reciprocal of the first 
term in the next, we have 


GLC. 


Now, to approximate a constant, we convert it to an integer 
fraction by letting a equal the constant stripped of its decimal 
point and carried to d digits, and letting b equal 10". As a 
pertinent example from Chapter 4, let us find rational approxima- 
tions to 10/1© = 1.154782.... Taking a = 1154782 and b = 
1000000, we proceed as follows: 


1154782 = 1(1000000) + 154782 
1000000 = 6(154782) + 71308 
154782 = 2(71308) + 12166 


| 
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71308 = 5(12166) + 10478 
12166 = 1(10478) + 1688 
10478 = 6(1688) + 350 
1688 = 4(350) + 288 
Stopping at this point, we have obtained the partial quotients 


a, = 1, a) = 6, a3 = 2, ay = 5, as = 1, ag = 6, and a7 = 4. 
These are used to determine the principal convergents A,: 


A; =a, = 1 
l 

A, = a, + — = 7/6 
a? 


Actually, at this point we can ease our work by realizing that 


for An = Palm 
Pn = 4nPp-1 = Pn-2 
An = 4nGn-1 * Gn-2 


and we can continue quite rapidly: 


peee tle 15 
D oe 
Dime G4] s 13 
58S at 7 ie 
" j20 +s 
97 (664 2753 
te —| ei: te = — 
5 ag MOm Sass AT 58, 


Now there are intermediate convergents as well between A,_> 
and A,, for each a, > 1. These are of the form: 
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Pn-2 


> ja 


Pn-2 + Pn-1 


Gn-2 + nel 


Pr + 2pm 
Qn-2 7 ler 


Pn-2 BP 4nPn-1 


Qn-2 7 An dn-1 


= A, 


Representing the intermediate convergents as (A, A,), 
we find 


2: (AAs) 1+7 8 
a= c — = — 
3 13 jee 7 
Pela e22. 37. 52 ar 
on 5: (A2A4) = = 


6+13 19°32’ 45'58 
as = 1: no intermediate convergents as a, < 1 


SZ +97 179 216 315 490 3G" 


71 + 84 155 239° 323 407’ 491 


ag = 6: (A4A6) = 


97 + 664 7/61 1425 2669 


o = 4: (ee eee 
q 84 +575 659’ 1234’ 1809 
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Here we find those principal convergents we would get if we 
allowed subtractions instead of additions in various steps of 
Euclid’s Algorithm. 

In the end, we arrive at two convergent sequences. The first one 
increases, approaching 10!/! from below, and consists of the odd 
principal convergents and their intermediate convergents: 


Ai; (A,A;3), A3; As; (AsA7), A? 
or, 


8 15 97 fet 25 2089 2753 
"7 13’ 84 659° 1234’ 1809 ’ 2384 
The second sequence decreases, approaching 10!/'© from above, 


and consists of the even principal convergents and their inter- 
mediate convergents: 


A), (A2A,), A4 (AyAg), Ag 
Of, 


7 22 By G2 Ge 179 276 373 470 367 GM 
6 19 32 45°58 71 155 239’ 323 407 491 | «575 


A computer provides, of course, a convenient means of tabulating 
these fractions, as it did for me. 

Now we have the intriguing task of choosing a fraction that 
provides good accuracy with a convenient numerator and denomi- 
nator for multiplication (we flip the fraction for division). Since 
a7, a4, and ag are relatively large, we are clued to look at A,, Aj, 
and As; for proximity to 1.154782 relative to their complexity. 
We discover 


A = | 
A; = 15/13 = 1.1538462... 
As = 97/84 = 1.1547619... 
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Perusing other choices reveals 52/45 = 104/90 as a possibility, 
but its accuracy is poor (1.1555555...). Therefore, perhaps our 
best pick, particularly in light of the factoring available for can- 
cellation when additional constants are multiplied or divided, is: 


10116 ~ 97/84 


100 
3 2 ew 


= 1.154762 compared to 1.154782... 


This rational approximation appears in Table 8. 
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Addition, 9-10 method with table for, 163-64 
Aitken, Alexander Craig, 5, 35, method without table for, 
77, 79 161-63 

Algebraic factors, 54-56 
Antilogarithms, 133ff 

Bemer’s method for, 137-42 

Brigg’s method for, 133-37 

definition of, 114, 113 
Arccosine, 156-62 

Borchardt’s Algorithm for, 


Bidder, George Parker, 117 

Borchardt’s Algorithm, 130-33, 
161, 164 

Briggs, Henry, 116, 133-34 


160-62 Chebyshev, P. L., 86 
method with table for, 157-58 | Composite numbers, 47. See also 
method without table for, Prime numbers. 
156-57 Computers 
precise method without table Bemer’s logarithm method on, 
for, 159-160 ey | 
Arcsine, 156-62 binary, 3, 38, 18, 132 
Borchardt’s Algorithm for, decimal, 3, 96, 137 
160-62 factoring methods on, 68 
method with table for, 157-58 hexadecimal test for, 45 
method without table for, precomputing constants with, 
166257 Lideh 
precise method without table Continued fractions, 172-78 
for, 159-60 definition of, 172 
Arctangent, 161-64 rational approximations using, 
Borchardet’s Algorithm for, 164 173-78 
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Cosine, 149-52 
method with table for, 151 
method without table for, 
149-51 
Cournot, A., 171 
Cross multiplication, 11 
Cube roots, particular methods 
for, 99-104. See also Roots. 
Cubing, 17-18. See also 
Multiplication. 


Degrees, conversion to radians, 
146, 172 
Divisibility, 39ff 
by numbers with multiple near 
power of ten, 41-42 
by odd numbers, 42-44 
by particular numbers, 40-42, 
45-46 
elevens test for, 40 
Euclid’s Algorithm for, 44 
factoring sieves using, 72-73 
nines test for, 39-40 
other bases for, 45—47 


Division, 19ff. See also Divisibility. 


by factors of numbers near 
round numbers, 28 

by numbers near round 
numbers, 24—28 

a cross (or Fourier), 28-33, 97 

repeating decimals and, 


19-24, 35 


Error checking, 39ff. See also 
Divisibility. 
Euclid’s Algorithm. See also 
Greatest common divisor. 
continued fractions and, 
173-78 
definition of, 35-38 
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divisibility tests with, 4445 
factoring with, 65 
Euler, factoring methods of, 
61-65. 


Factorization, 47ff 
Euler’s method of, 61-65 
Fermat’s method of, 48—54, 
68-73 
general Euler’s method of, 65 
Legendre’s method of, 54-56 
Triangular number method of, 
57-61, 65-66, 73 
Vaes’ method of, 56-57 
Fermat, Pierre de, 48, 54, 61 
Fifth roots, particular methods 
for, 105-108. See also Roots. 
Finger counting, 12 


Fourier, Joseph, 28, 97 


Greatest common divisor, 35ff 
Euclid’s Algorithm for, 35-36 
halving method for, 38 
least-remainder method for, 36 
modified Euclid’s Algorithm 

fot; 37-365 


Halley, Edmond, 83, 107, 123 

Hexadecimal error checking, 
45-46 

Hofstadter, Douglas R., 4-5 

Hyperbolic trigonometric 
functions, methods for, 165 


Inverse logarithms. See Anti- 
logarithms. 
Inverse trigonometric functions, 


156ff. See also Arccosine; 


Arcsine; Arctangent. 
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Keyser, C. J., 167 
Knuth, Donald E., 18, 48 
Kogbetliantz, E. G., 113 


Lagrange interpolation, 129-30 
Legendre, A. M., factoring 
method of, 54-56 
Lightning calculators 
methods of, 11, 18, 20-21, 24, 
35, 44, 62, 77-79, 117-18 
names of, 4-5, 117 
performance of, 3-5 
Logarithms, 113ff 
Bidder’s method for, 117-18 
Borchardt’s Algorithm for, 
130-33 
definition of, 114-15 
factoring for, 40 
Halley’s method for, 123-26 
Lagrange’s formula applied to, 
129-30 
neighboring value methods for, 
121-23 
Newton’s method for, 126-28 
prime number, 119-21 
root extraction using, 142-144 
series approximation for, 
116-17 
Lowell, James Russell, 1 
Lowest common multiple, 39 


Melding, 6-8 
Mental calculators. See Lightning 
calculators. 
Modular arithmetic 
identities in, 35-36 
sieves using, 50-56, 59-61, 64, 
68-73 
Multiplication, 10-18. See also 
Cubing; Squaring. 
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cross, 11 

distributive law method of, 11 
large number, 18, 47 
particular simplifications in, 


12-15 


Napier, John, 113, 116 
Newton, Isaac, 87, 101, 126 
Notation 

mod, 26, 36 

pi, 129 

sigma, 29, 129 

vertical bar, 6-8, 10, 89 


Ozanam, triangle relation of, 152, 
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Prime numbers. See also 
Factorization. 
definition of, 19, 48 
fundamental theorem of 
arithmetic and, 47-48 
logarithms of, 118-121 rough 
test for, 66—68 
sums of squares of, 61-62 
Primitives, 9ff 
definition of, 9 


Quadratic residues, 50 


Radians, conversion to degrees, 
146, 172 

Rational approximations, 171 ft 

Reciprocals. See also Division. 
approximations to, 33-34 
calculation of, 19-24 
square roots of, 97-99 

Roots, 77ff. See also Cube roots; 
Fifth roots; Square roots. 
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Chebyshev’s correction for, 
86-87, 98-99, 100-102, 
105 

general algorithm for square, 
87-97 

Halley’s method for, 83-85, 98, 
100, 105 

integer, 77-79 

Newton—Raphson method for, 
80-83, 97-99, 103-105 

particular methods for cube, 
99-104 

particular methods for fifth and 
higher, 105-108 

particular methods for square, 
79-87, 98 

reciprocal square, 97-99 

use of logarithms for, 115, 
142-44 

Russell, Bertrand, 145 


Sine, 145-52 
method with table for, 146-48 
method without table for, 
148-149 
Smith, Steven B., 3-4, 44, 79 
Square roots, 77-99. See also 
Roots. 
approximation for reciprocal, 


97-99 
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general algorithm for, 87-97 
particular methods for, 79-87 
Wallis’ extraction of, 109 
Squares 
endings of, 49 
endings of differences of, 69 
sieves for, 50-52, 69-70 
Squaring, 15-17. See also 
Multiplication. 
Subtraction, 9-10 
Sylvester, J. J., 167 


Tangent, 152-156 
general formula for, 153—54 
precise formula for, 154-56 
sine/cosine method for, 152 
Triangular numbers 
definition of, 57-58 
endings of, 59 
factoring methods using, 
57-61, 65-66, 73 
sieves for, 59-61 
Trigonometric functions, 145ff. 
See also Cosine; Sine; 
Tangent. 


Vaes, factoring method of, 56-57 


Wallis, John, 5, 109 
Whitehead, Alfred North, 9 
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“Fascinating . . . no author has gone as far as 
Doerfler in covering methods of mental 
calculation beyond simple arithmetic. Even 
if you have no interest in competing with 
computers you'll learn a great deal about 
number theory and the art of efficient 


computer programming.” —Martin Gardner, 
Author of Mathematical Carnival 
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